Tumors from 24 adult tissue types separate into two distinct clusters based on expression of epifactor genes

We investigated whether tumors from each of the 24 adult cancer types in the TCGA repository (Fig.1a and Supplementary Data1)16 would separate into well-defined subgroups based on the expression patterns of 720 epifactor genes from the Epifactors database23 (Supplementary Data2). These epifactors encode proteins involved in the addition, removal, and recognition of DNA methylation and histone marks, and chromatin remodeling (Fig.1b, top panel,Supplementary Data2). The majority of the epifactor genes (556 out of 720) are not known to be genetically altered in cancer tissues (Fig.1b, bottom panel and Supplementary Data2)24,25. We clustered the patient tumors from each cancer type using the non-negative matrix factorization (NMF) algorithmbased on the epifactor genes with the most variable expression among the patient tumors (Fig.1c and Supplementary Data 1)26,27. With NMF clustering, a reduced representation of the gene expression data is generated that delineates a subset of genes that are important for separating the samples into clusters. For each of the 24 cancer types evaluated independently, separating the tumors into two clusters resulted in the best solution based on three measures of cluster validation (Supplementary Data1, Supplementary Fig.1, and Supplementary Data3). The two clusters for each cancer type were characterized by a set of signature top NMF genes with distinct expression patterns for the tumors in the two clusters (Supplementary Data4). As an example, for breast cancer, two distinct tumor clusters were observed (Fig.1d and Supplementary Fig.2a), and the PAM50 breast cancer subtypes28,29 were non-randomly distributed between the two BRCA epifactor expression-based clusters (Supplementary Fig.2b), consistent with a previous study showing different epigenetic characteristics for the PAM50 subtypes30.

a Tissue of origin for the 24 adult cancer types from TCGA included in the clustering analysis based on epifactor expression. Tissue locations are labeled with TCGA abbreviations. Sex-specific tissue locations are shown in purple for female (left panel) and blue for male (center panel). The full names for each cancer type are provided (right panel). b Functional categories for the 720 epifactor genes included in this study (top panel). This list of epifactors genes was obtained from the manually curated Epifactors database generated by Medvedeva et al.23. The bottom panel shows the overlap of these epifactor genes with the genetically altered cancer genes cataloged in either the COSMIC24 or OncokB25 databases. c The NMF-based clustering26,27 analysis workflow is provided. Raw RNA-seq counts for all of the genes in each patients tumor for a specific cancer type were normalized using the DESeq267 R package. The most variable epigenetic genes (var_epi) were selected based on a cancer type-specific standard deviation cutoff. This dimensionally reduced counts matrix (patient x var_epi) was used as an input to the NMF R program27. d PCA plot showing the two clusters (red and cyan) of the BRCA patient tumors (depicted as dots) as determined by the NMF method. The variances explained by principal components (PC) 1 (x axis) and PC2 (y axis) are plotted. e Heatmap showing the fraction of the top NMF genes for the cancer type in the corresponding column that overlaps with the top NMF genes for the cancer type in the corresponding row. Darker colors indicate a higher fraction of overlap. The rows and columns are hierarchically clustered. f Heatmap describing the most frequent top NMF genes (rows, genes ranked in decreasing order of frequency) across the 24 cancer types (columns). a and b were created using Biorender. Supporting information for this figure can be found in Supplementary Figs.1 and 2, and Supplementary Data14.

The number of top NMF genes across the 24 cancer types ranged from 76 genes for LGG to 9 genes for CRC, with a median of 43 genes (Supplementary Fig.2c). A pan-cancer map based on the expression patterns of the top NMF genes from all tumor types showed that the tumors group largely based on their tissues of origin, and, to some extent, tissue proximity (Supplementary Fig.2d). For example, KIRC and KIRP, two types of kidney cancer, were found near each other, and LGG and GBM, two types of brain cancer, were also adjacent to each other in this low-dimensional representation.

There was a high overlap among the top NMF genes in the ACC (carcinoma of the adrenal glands that sit atop each kidney), BRCA, LIHC (liver), LUAD (lung), STAD (stomach), LGG (brain), UCEC (uterus), and SARC (soft tissues and bone) cancer types (Fig.1e). A strong overlap was also observed among the top NMF genes in the KIRC (kidney), KIRP (kidney) and PAAD (pancreas) tumors (Fig.1e). SATB1 was the gene most frequently represented as a top NMF gene and was a signature gene for 12 cancer types (Fig.1f). SATB1 mediates chromatin organization by acting as a landing platform for chromatin remodeling proteins31. When the top NMF genes were assigned to one of 19 protein complex groups23, epifactors belonging to histone acetyltransferase (HAT) complexes were significantly enriched (P <0.05)among the top NMF genes for seven cancer types (LIHC, GBM, UCEC, BRCA, LUAD, KIRP, and PAAD).

To determine whether patients in the two clusters developed from expression levels of epifactors differ with regard to their clinical outcomes, we compared the progression-free interval (PFI), disease-specific survival (DSS), and overall survival of the patients in the two clusters for each cancer type. Cox regression was used to adjust for the effects of age and sex of the patients in each cluster, unless otherwise mentioned. The clusters from 10 out of 24 cancer types (ACC, CRC, KIRC, KIRP, LGG, LIHC, LUAD, PRAD, STAD, and UCEC) had significant differences in clinical outcome (P<0.05) for at least one of the three metrics (PFI, DSS, and overall survival) (Fig.2a).

a Heatmap showing the significance (P value from multivariate Cox regression analysis; adjusted for age and sex) of the difference in the clinical outcome (PFI, DSS, and overall survival) between the two epifactor expression-based tumor clusters for each of the 24 cancer types. The grey color indicates that the difference in clinical outcome between the two clusters is not significant. bf KaplanMeier plots comparing the progression-free intervals of the two NMF-derived clusters for the five-cancer group that show significant differences in clinical outcome for the three metrics PFI, DSS, and overall survival. Significance was determined with the log-rank MantelCox test. The cluster with poor outcome is designated (in superscript) poor, while the cluster with better outcome is designated better. The number of patients (n) in each cluster is shown. b ACCpoor n=40, ACCbetter n=31. c KIRCpoor n=61, KIRCbetter n=108. d LGGpoor n=107, LGGbetter n=146, e LIHCpoor n=70, LIHCbetter n=72, f LUADpoor n=20, LUADbetter n=49. g Heatmap showing the significance (P value; two-tailed Fishers exact test) of the difference in clinical metrics (pathologic M, pathologic T, pathologic N, stage, and grade) between the epifactor expression-derived clusters for the five-cancer group. hl Barplots of the clinical characteristics for instances (shown in g) in which the two clusters significantly differ. m, n Composition of epifactor expression-derived clusters for ACC (m) and LGG (n) with regard to established TCGA subtypes. o Classification of the epifactor expression-derived clusters for the five-cancer group based on established immunologic subtypes from ref. 33 with the following prognostic order (worst to best): C4~C6>C2~C1>C3~C5. Data for the clinical metrics were obtained from cBioPortal for Cancer Genomics60. Supporting information for this figure can be found in Supplementary Figs.38 and Supplementary Data5.

For the cancer types in the five-cancer group (ACC, KIRC, LGG, LIHC, and LUAD), the two clusters (Supplementary Figs.3 and 4) significantly differed in clinical outcome for all three metrics (Fig.2bf and Supplementary Fig.5). Consistent with the differences in outcome, the poor outcome ACCpoor cluster was composed of tumors with higher cancer stage (TNM stages 3 and 4), larger size (T3 and T4), and greater likelihood of lymph node spread (N1), compared to tumors in the better outcome ACCbetter cluster (Fig.2gj). The poor outcome LGGpoor and LIHCpoor tumors also included a significantly higher fraction of patients with grade 3 tumors than the tumors in the LGGbetter and LIHCbetter clusters, respectively (Fig.2g, k, l). When the clinical outcome differences between the NMF clusters were adjusted for stage and grade, all three metrics were still significant, except that for LGG and LIHC, significance was observed for 2 out of 3 metrics (Supplementary Fig6a). The distributions of the patients races and ethnicities did not differ between the clusters (Supplementary Fig6bg and Supplementary Data1) except for LIHCpoor, which had a higher fraction (~twofold) of Asian patients than LIHCbetter.

The prognostic efficacy of epifactor expression-based clusters was better than grade or epithelial-to-mesenchymal transition (EMT) for the five-cancer group (Supplementary Fig.7ad). Tumor grade was effective in predicting the outcome for just 1 or 2 cancer types, out of the five, across the three outcome metrics, while EMT could predict outcome for 24 cancer types across the three outcome metrics (Supplementary Fig.7ad). The two tumor clusters were not significantly different for EMT for ACC, LGG, LIHC, or LUAD (Supplementary Fig.7e). For KIRC, we did observe a significant difference (P=0.0015, two-tailed MannWhitney test) in the EMT scores between the two clusters, however, the tumors in the poor outcome cluster had a lower EMT score (median = 0.0750) compared to the tumors in the better outcome cluster (median = 0.17), the opposite of the expectation that a more mesenchymal phenotype would be associated with a worse prognosis32.

We asked whether the clinical differences between the clusters might reflect differences in the purities of the tumors in the two clusters. There was a significant difference in purity only for ACC and LGG (P<0.05, two-tailed MannWhitney test) (Supplementary Fig.8a), and all three clinical outcomes (PFI, DSS, and overall survival) were still significantly different for both ACC (P=0.002, P=0.006, and P=0.001; Cox regression) and LGG (P=0.008, P=0.007, and P=0.02) after adjusting for tumor purity. The differences in stromal fraction were significantly different for ACC and LGG clusters (Supplementary Fig.8b), but the clinical outcome differences between the NMF clusters for ACC (P=0.0003, P=0.0006, and P=0.0003; Cox regression) and LGG (P=0.0006, P=0.001, and P=0.003) were still significant after adjusting for stromal fraction. The ACC and KIRC clusters had different levels of immune infiltration (Supplementary Fig.8c), but the clinical outcome differences between the NMF clusters for ACC (P=0.007, P=0.015, and P=0.002; Cox regression) and KIRC (P=0.003, P=0.005, and P=0.0007) were still significant after adjusting for leukocyte infiltration.

We compared the epifactor expression-derived clusters with established TCGA subtypes for the five-cancer group33. None of the clusters were composed of only a single TCGA-defined tumor subtype. But, for each of the clusters, there was at least one TCGA subtype that was overrepresented (Fig.2m, n, Supplementary Fig.8df and Supplementary Data5). For example, for ACC and LGG (Fig.2m, n), there was a greater representation of some DNA methylation-based TCGA subtype(s) in the poor outcome cluster compared with the better outcome cluster, and vice versa. These results indicate that previously reported TCGA subtypes may have epigenetic features that contribute to their distinctive characteristics. Our epifactor expression-derived clusters also contained tumors with significantly different compositions of immunologic subtypes (see Methods)33 (Fig.2o and Supplementary Data5). Epifactor expression-derived clusters with poor clinical outcomes were enriched in immunological subtypes associated with poor prognosis (such as C4) and/or depleted of the subtypes associated with better outcome (such as C3 and C5), consistent with epifactor expression in cancer cells affecting the immune response to the tumor.

We performed a detailed analysis to determine the significant differences (adjusted P value<0.05, BenjaminiHochberg method) in the frequencies (fraction of affected patient tumors) of mutations and copy number alternations (CNAs) between the two clusters for these five cancer types (Supplementary Data4). For ACC, there were no significant differences in the mutation or CNA frequencies between the clusters for any gene (epifactor or non-epifactor). For KIRC, the two clusters were different in terms of CNA frequencies for six epifactor and 431 non-epifactor genes. None of these six epifactors (NPM1, UIMC1, NSD1, HDAC3, DND1, and TAF7) were assigned as cluster-defining top NMF genes. No differences in mutational frequencies between the clusters were observed for any gene. For LGG, only three non-epifactor and zero epifactor genes had significant mutational frequency differences, while three epifactor and 58 non-epifactor genes had CNA frequency differences, between the two clusters. The three epifactors (PRMT8, CHD4, and ING4) with CNA frequency differences were not part of the cluster-defining top NMF gene group. For LIHC, none of the genes had a difference in mutational frequencies between the two clusters. Twenty epifactor genes including one top NMF gene (TONSL), and 790 non-epifactor genes had significant differences in CNA frequencies between the two clusters. The CNA in TONSL affected <16% of patient tumors in both the clusters suggesting that this CNA is unlikely to have a major effect on the observed differences in patient outcome. For LUAD, only TP53, an epifactor gene, displayed differences in mutational frequencies between the two clusters (P=0.006; 68% in LUADpoor vs. 21% in LUADbetter, two-tailed Fishers exact test and adjusted for multiple hypothesis correction). None of the genes (epifactor or non-epifactor) showed differences in CNA frequencies between the two clusters. After adjusting for TP53 mutations, all three clinical outcomes were still significantly different for the two LUAD clusters. These results suggest that the differences in expression levels of signature epifactor genes for the poor and better outcome clusters were unlikely to exclusively reflect mutations or CNAs.

We performed weighted correlation network analysis (WGCNA)34 to identify the gene ontology (GO) terms35 associated with gene groups (modules) with similar patterns of expression as the top NMF epifactor genes of poor outcome or better outcome clusters (Fig.3a, b, Supplementary Fig.9ac and Supplementary Data6). The GO terms for the modules related to poor outcome clusters were enriched for cell cycle genes (dark orange module, ACC; midnight blue module, LGG; grey60 module, LIHC; and green module, LUAD) and developmental genes (turquoise module, LIHC), indicating that differences in proliferation rate or stem-like features36 may contribute to the clinical differences observed between the clusters. The protein-protein interaction (PPI) networks formed from the top NMF epifactors were significantly enriched (P<0.05) compared to background for all the five cancer types (Supplementary Data6, Fig.3c, d and Supplementary Fig.9df). The top NMF epifactors belonging to the cell cycle-related modules formed tight, well-connected PPI networks indicating a possible coordinated mechanism of action37.

a, b GO terms (bar labels on the right) associated with gene modules containing top NMF genes for ACC (a) and LIHC (b). Modules were generated using the WGCNA analysis tool34 applied to co-expressed genes. Only modules containing at least five top NMF genes were considered. Modules in which top NMF genes are associated with poor outcome are shown in blue, and modules for better outcome are shown in purple. Adjusted P values (Padj) were obtained by applying BenjaminiHochberg multiple test correction to the unadjusted P values in a module. Only the top representative GO terms related to biological process or molecular function were considered for each gene module. c, d PPI networks were generated for the encoded proteins of the top NMF genes for the ACC (c) and LIHC (d) cancer types. The top NMF genes in the networks (nodes; shown as circles) are colored based on the modules in which they reside. Top NMF genes that were not assigned to a module with 5 or more top NMF genes were depicted as white circles (no color fill). The thickness of a line (edge) connecting two top NMF genes (nodes) indicates the confidence level of the protein-protein interaction prediction between those two top NMF genes. Supporting information for this figure can be found in Supplementary Fig.9 and Supplementary Data6 and 7.

PCA plots based on array-based DNA methylation levels for the five cancer types13 revealed that DNA methylation captures some of the differences between the clusters developed based on epifactor gene expression, but that DNA methylation alone provides significantly less separation between the two tumor clusters than can be achieved by analyzing data from all epifactor genes (Supplementary Fig.10).

DNA methylation factors tend to be expressed at higher levels in tumors with poor outcome (red color in the heatmapinSupplementary Fig.11a). This is true for the epifactors that are directly involved in de novo DNA methylation (DNMT3A and DNMT3B) or in the maintenance of DNA methylation (DNMT1 and UHRF1)13. For each tumor type, we determined the number of hypermethylated and hypomethylated loci in the poor outcome cluster compared to the better outcome cluster (Supplementary Fig.11b and Supplementary Data7). The pattern of differential methylation between the clusters varied across the five cancer types with more hypermethylation events in the tumors in poor outcome clusters for ACC and LIHC, while the reverse was true for LGG. We determined for all hypermethylated and hypomethylated sites linked to genic regions, whether the closest gene was upregulated or downregulated in the poor outcome cluster (Supplementary Fig.11cg). The relationship between DNA methylation state and gene expression levels was significant (Fishers exact test) for all cancer types except KIRC, suggesting an impact of DNA methylation levels on downstream gene regulation.

To complement our clustering analysis based on patterns detected among all of the epifactors, we performed a systematic analysis of the prognostic value of the expression levels of each of the variable epifactor genes (see Methods) considered individually across the 24 cancer types (Supplementary Data8). The fraction of prognostic epifactors (out of the total number of variable epifactors) varied across the cancer types and ranged from 77% for KIRC to 0.4% for TGCT (Fig.4a), with a median of 21%. The prognostic direction of a gene was not always the same across the cancer types (Supplementary Data8) and the expression levels of these prognostic genes among the tumors were not consistently associated with mutations or CNAs that could explain the expression level differences associated with patient outcome (Supplementary Data8). Among the 24 cancer types, the fractions of prognostic epifactors and non-epifactors were highly correlated (Supplementary Fig.12a, b), but on average, the fraction of prognostic genes was higher for epifactor genes than for non-epifactor genes (P=0.015, Wilcoxon matched-pairs signed rank test) (Supplementary Fig.12c). The fraction of variable epifactors that are prognostic among tumor types had a weak negative correlation that did not reach statistical significance (P>0.05) with either the total number of mutations (Supplementary Fig.13a) or the total number of copy number alterations (CNAs) (Supplementary Fig.13b).

a Number of prognostic epifactor genes for each of the 24 cancer types. Prognostic genes were identified based on a significant difference in PFI outcome between patient tumors with high and low expression levels of the gene. p values were adjusted for age and sex of patients and for multiple hypothesis testing (BenjaminiHochberg method). b Circos plot showing the epifactor genes that are most frequently prognostic across cancer types. The lines connect the genes and the cancer types in which they were determined to be prognostic. Red lines indicate the five-cancer group (ACC, KIRC, LGG, LIHC, and LUAD) and blue lines indicate other cancer types. c Forest plot showing the hazard ratio and 95% confidence interval (CI) for the most significant prognostic gene (gene symbols in parentheses) for each of the 24 cancer types. Hazard ratios lower than one indicate that higher expression of the gene is associated with poorer PFI. d Heatmap for the enrichment of protein complexes among the prognostic genes for the 24 cancer types. White rectangles indicate no significant enrichment. Significance was determined using permutation tests. e Barplot showing the fraction of top NMF genes that are prognostic for the 24 cancer types. f Heatmap indicating the prognostic status of the most frequent top NMF genes across the five cancer types. g KaplanMeier plots for the most significant prognostic SWI/SNF genes for ACC and LIHC. Significance was determined with a log-rank MantelCox test and the number of patient tumors (n) in each group are provided. ACCSMARCD1,high n=20, ACCSMARCD1,low n=51, LIHCARID1A,high n=19, LIHCARID1A,low n=123, LIHCCHAF1B,high n=60, LIHCCHAF1B,low n=82. h Bar plots indicating the effect of meta-PCNA correction on the number of cancer types for which a gene is prognostic for the top NMF genes that were included in WGCNA-derived gene modules for ACC (left and center) and LIHC (right). i Prognostic status for genes that remain prognostic after the meta-PCNA correction among the cancer types. Supporting information for this figure can be found in Supplementary Data8.

The top ten most frequent prognostic epifactor genes across the cancer types (Fig.4b) were involved in chromatin remodeling (DPF1 and TOP2A) and in depositing and reading histone modifications including histone phosphorylation, methylation, and deubiquitination (AURKA, BUB1, CDK1, CHEK1, GSG2, MSH6, SMYD2, and USP49). Out of these, AURKA, TOP2A, CDK1, and BUB1 were also included in the list of most frequent top NMF genes across the 24 cancer types (Fig.1f). For the most significant prognostic gene for each of the 24 cancer types, high expression of the prognostic gene was associated with poor outcome for 15 cancer types (hazard ratio <1), while high expression of the prognostic gene was associated with better outcome for nine cancer types (hazard ratio >1) (Fig.4c). Genes associated with the HAT or chromatin remodeling (SWI/SNF or ISWI) complexes were significantly overrepresented among the prognostic genes (P<0.05) in 11 or more cancer types (Fig.4d).

For each epifactor, we determined the number of cancer types in which high expression (N) or low expression (M) of that epifactor was associated with poor outcome. Positive prognostic residuals (N-M>0) were more frequent than negative prognostic residuals, indicating that high expression of epifactors was more often associated with poor outcome for the epifactors overall (all groups in Supplementary Fig.14a, b) and for subsets of epifactors associated with DNA modification (n=25), histone modifications (n=487), and chromatin remodeling (n=124) (Supplementary Fig.14a, b). We observed a similar trend among epifactors that are histone writers (n=147), erasers (n=57), and readers (n=80) (Supplementary Fig.14c, d). When we further divided the histone writers based upon the specific histone mark they deposit, the writers that catalyze histone acetylation (n=33), but not those that catalyze methylation (n=47), phosphorylation (n=36), or ubiquitination (n=28), had more negative prognostic residuals than positive residuals, indicating that high expression of histone acetylases is associated with a more favorable prognosis (Supplementary Fig.14e, f).

A higher fraction of the top NMF genes for the five-cancer group were prognostic for outcome as compared with other cancer types (Fig.4e, f). Further, the prognostic direction of these frequent top NMF genes was consistent across the five-cancer group, with the exception of PPARGC1A (Fig.4f). The frequent top NMF genes ASF1B, ATAD2, BUB1, CDK1, CHAF1A, HJURP, PBK, and TOP2A (Fig.4f) that were signature genes for the poor outcome cluster in ACC, LGG, LIHC, and LUAD (Supplementary Fig.4 and Supplementary Data4) were also significantly associated with poor outcome when expressed at high levels for those same cancer types (Fig.4f). The prognostic genes for ACC and LIHC had a shared enrichment for SWI/SNF chromatin remodeler genes (Fig.4d) with SMARCD1 in ACC, and ARID1A and CHAF1B in LIHC, being the most significantly prognostic SWI/SNF genes (P<0.0001) in these two cancer types (Fig.4g). The top NMF epifactors were more likely to be predictive of outcome than SWI/SNF epifactors, overall (Supplementary Fig.15a, b).

In their investigation of genes that predict breast cancer outcome, Venet et al. found that the prognostic value of the majority of signature genes was eliminated when they adjusted for the expression levels of a meta-PCNA signature that removed the confounding effects of cell proliferation38,39. After adjusting for the meta-PCNA signature, in addition to age and sex, we found that the expression levels of some prognostic epifactor genes were no longer associated with clinical outcome (Supplementary Data8). Analysis of the instances in which an epifactor gene was proliferation-independent (prognostic even after meta-PCNA correction) or proliferation-dependent (not prognostic after meta-PCNA correction) revealed that the top NMF genes belonging to a co-expression module (Fig.3a, b and Supplementary Fig.9ac) with a highly enriched cell cycle GO term (such as the dark orange module of ACC) were more affected by meta-PCNA correction compared to top NMF genes in modules highly enriched for other GO terms such as autophagy (the blue module of ACC) or development (the turquoise module of LIHC) (Fig.4h). The frequently prognostic epifactor genes that were also unaffected by the meta-PCNA correction included those involved in histone modifications (KDM4B, KAT6A, MBTD1, MTA1, and PHF1), histone binding (BRD3, MBTD1, and PHF1), and chromatin organization and remodeling (MTA1) (Fig.4i).

We asked whether pan-cancer epigenetic features can be used to develop a predictor for patient outcome for the five-cancer group. To achieve this, we used the Cox-nnet artificial neural network (ANN) framework by Ching et al.40. The Cox-nnet model consists of an input layer, a hidden layer with 143 nodes, and a final Cox-regression layer that outputs the prognostic index (PI), equivalent to the log hazards ratio (Fig.5a). Patients from the combined cohort of the five cancer types were randomly split (80:20) into training and test sets. For the model trained on the epifactor expression data, age and sex of the patients in the 5-cancer group, the clinical outcomes (PFI) for the high PI and low PI groups of the test set were significantly different (P<0.0002) (Fig.5b), indicating that the trained model was able to successfully predict the likely clinical outcome for patients that were not included in the training set. As the top NMF epifactor genes of KIRC showed less overlap with the remaining four cancer types (Figs.1e and 4f), we also trained a Cox-nnet model based only on the other four cancer types (ACC, LGG, LIHC, and LUAD). With this 4-cancer-type model, the log-rank P value for the test set was highly significant (P<0.0001) (Fig.5c). The model trained only on KIRC did not result in groups with a significant difference in outcome (P=0.19) (Fig.5d).

a A Cox-nnet model40 was used as a framework for predicting patient outcomes. The patient cohort was randomly split (80:20) into training and test sets. The model was trained on input features consisting of the expression values of the 720 epifactor genes, and the age and sex of the patients in the training set. The model consisted of an input layer that accepts the input features and is fully connected to a hidden layer. The output of the nodes of the hidden layer was fed to a cox-regression layer. The final output of the model was the log hazard ratios of the patients (prognostic index, PI). To evaluate the performance of the model, the test set patients were divided into high PI and low PI groups based on the median PI of the patients. The clinical outcomes between these two groups were compared using the log-rank MantelCox test (KaplanMeier method). Created using Biorender. bd KaplanMeier plots evaluating the performance of the model. b Results when the model was trained and tested on patients from the 5-cancer group (ACC, KIRC, LGG, LIHC, and LUAD). High PI n=71, Low PI n=70. c Results when the model was trained and tested on four cancer types (ACC, LGG, LIHC, and LUAD). High PI n=55, Low PI n=56. d Results for a model trained and tested on only KIRC. High PI n=17, Low PI n=17. e Prognostic status of the top 20 input features (left panel) ranked on the basis of their importance in the Cox-nnet machine learning (ML) model for the five cancer types is shown. A heatmap indicating which of the top 20 features from the left panel are also top NMF genes across the five cancer types is shown on the right. Only the features that are a top NMF gene for at least one cancer type are shown. f Same as (e), but for the four cancer type model. KIRC was not included in the Cox-nnet model, but is included in these heatmaps for comparison. Supporting information for this figure can be found in Supplementary Data9.

Most of the top 20 important features for clinical outcome from the pan-cancer model (Supplementary Data9) were individually prognostic (P<0.05; Supplementary Data8) with higher expression of these features associated with poor outcome (Fig.5e, f, left panels). About half of these important features (10 out of 20 for the 5-cancer model and 11 out of 20 for the 4-cancer model) were top NMF genes in at least one of the 5 cancer types (Fig.5e, f, right panels).

To further test and validate our findings on the prognostic role of epifactors, we used independent, publicly available datasets for KIRC, LGG, and LUAD (Supplementary Data9). For each cancer type, we assigned the tumors in this validation cohort to either poor outcome or better outcome groups (Supplementary Fig.16ag and Supplementary Data9) based on the expression pattern of the top NMF epifactor markers that we determined based on the original datasets (Supplementary Fig.4). In the case of KIRC and LUAD, we observed significant clinical differences (P<0.05; Cox regression, adjusted for age and sex) between the two groups of tumors (Supplementary Fig.16bg), while for LGG, the difference was nearly significant (P=0.071) (Supplementary Fig.16a). There was also a significant overlap (P<0.05, based on the hypergeometric distribution) between the epifactors that were individually prognostic for the validation and primary datasets (Supplementary Fig.16h and Supplementary Data9) for KIRC (P=0.019), LGG (P=0.0003), and LUAD (P=0.014). These results demonstrate that expression levels of these epifactors, together or individually, have a robust capacity to classify tumors based on clinical outcomes.

Mutation frequencies are estimated to be 14 times lower in pediatric than adult cancers41. In one detailed genomic study, for 10% of pediatric tumors, no underlying, cancer-promoting mutation or structural copy number variant could be identified41. From this perspective, pediatric tumors have the potential to be more epigenetically driven than adult tumors. To compare our findingson epifactors in adult tumors with pediatric cancers, we obtained genomic and clinical data for pediatric tumors from four high-risk cancer types (neuroblastoma (NBL), osteosarcoma (OS), acute myleoid leukemia (AML), and Wilms tumor (WT)) in the Therapeutically Applicable Research to Generate Effective Treatments (TARGET) datasets (Fig.6a). These four pediatric cancer types are difficult to treat and originate in different tissues and cells within the body: immature nerve cells of the sympathetic nervous system (NBL); bone (OS); immature white blood cells of the bone marrow (AML); and kidney (WT).

a Schematic depicting the four high-risk and hard-to-treat pediatric cancer types (NBL, OS, AML, and WT) from the TARGET program in this study. These cancer types originate in brain (NBL), bone (OS), immature white blood cells (AML), and kidney (WT) in children and adolescents. These pediatric tumors were compared with 24 adult cancer types (depicted in Fig.1a) from TCGA. b Two epigenetic expression-based clusters showed significantly different survival outcomes for NBL and OS (KaplanMeier survival plots). NBLpoor n=68, NBLbetter n=37. OSpoor n=40, OSbetter n=28. c Heatmap showing the fraction of cluster-defining top NMF genes for four pediatric cancer types (columns) that overlap with the top NMF genes of the 24 adult cancer types (rows). d Heatmap showing status of the most frequent pediatric top NMF genes (rows) as top NMF genes for the different pediatric and adult cancer types. e Enrichment of different protein complexes in the top NMF genes of the four pediatric cancer types. Only the significantly enriched complexes (purple) are shown. Significance was calculated using the two-tailed Fishers exact test. f Barplot indicating the GO terms for three different gene modules for NBL obtained from the WGCNA analysis. These modules included at least 11 top NMF genes that were either all upregulated in the poor outcome or better outcome clusters obtained from the NMF algorithm. g PPI network generated using the top NMF genes of NBL. The genes are color-coded based on their associated WGCNA modules (shown in f). Thicker edges (connecting lines between the genes (nodes)) indicate higher degree of PPI between the two genes connected by the edge. h KaplanMeier survival plots for the most prognostic gene for NBL (RUVBL2) and OS (PRDM12). NBLRUVBL2,high n=16, NBLRUVBL2,low n=89. OSPRDM12,high n=13, OSPRDM12,low n=55. i Heatmap indicating the prognostic value of the most frequent pediatric prognostic genes (rows) across the pediatric and adult cancer types (columns). White indicates the gene is not prognostic; green indicates that low expression is associated with a poor outcome; pink indicates that high expression is associated with poor outcome. For (b, h), the P values from the log-rank MantelCox and the number of patients in each cluster (n) are indicated. Supporting information for this figure can be found in Supplementary Fig.7 and Supplementary Data10 and 11.

NMF clustering based on expression levels of epifactors in the pediatric patient tumors resulted in two clusters for each of the four pediatric cancer types (Supplementary Data10). The two clusters were significantly different in overall survival (corrected for age and sex, Cox regression) for NBL (P=0.024) and OS (P=0.032) (Fig.6b and Supplementary Fig.17ad), but not for AML (P=0.092) and WT (P=0.693). The top NMF genes for the four pediatric cancer types (81 top NMF genes for NBL; 34 for OS; 27 for AML; and 31 for WT) (Supplementary Data10) overlapped to different degrees (ranging from 0 to 44%) with the top NMF genes of the 24 adult cancer types from TCGA (Fig.6c). Out of the 21 genes that were a top NMF gene in at least two pediatric cancer types (Fig.6d), seven genes (ASF1B, AURKB, SMARCD3, TONSL, UBE2T, ZBTB7C, and ZNHIT1) were shared with the 27 most frequent top NMF genes in adult cancers (Fig.1f). Across the 24 adult cancer types, LIHC and OV had the most overlap of 8 genes between their top NMF genes and the frequent pediatric top NMF genes (Fig.6d, right heatmap). The top NMF genes of the pediatric cancer types were enriched for genes related to SWI/SNF (WT and OS), HMT (WT), MLL (WT), and HAT (NBL) protein complexes23 (Fig.6e). Similar to our findings for the adult cancer types, the signature top NMF genes for the poor outcome clusters of both NBL and OS (Fig.6b) correlated in expression with cell cycle genes in the green modules of NBL (Fig.6f) and OS (Supplementary Fig.17e). Also similar to our findings in adult tumors, the top NMF genes included in the green module for NBL formed a well-connected PPI interaction network (Fig.6g).

Out of 720 epifactors, 51 genes were prognostic for overall survival in NBL, 98 genes in OS, and 97 genes in AML (Supplementary Data11). None of the epifactor genes were prognostic for overall survival in WT. KaplanMeier plots for the most significantly prognostic genes for NBL (RUVBL2) and OS (PRDM12) are shown in Fig.6h. The prognostic value did not change after the meta-PCNA correction for any of the pediatric prognostic genes (Supplementary Data11). Twenty-four epifactor genes were prognostic in at least 2 of 3 pediatric cancer types (Fig.6i, left heatmap). The prognostic value (and direction) of these 24 epifactor genes varied across the adult cancer types (Fig.6i, right heatmap) with the most overlap observed for ACC, KIRC, and LGG, three members of the 5-cancer group.

Given our observation that tumors were associated with one epifactor expression-based cluster or another, we asked whether the individual cells within a tumor would display a gene expression profile related to one of the two clusters. Expression-based composite scores derived from the signature genes for each of the two epifactor expression-based clusters (Supplementary Fig.3c) were mapped to each cancer cell in a two-patient, single-cell RNA-seq dataset for LGG42 (see Methods). Individual cancer cells were assigned to one of the four different groups: LGGpoor, LGGbetter, LGGpoor+LGGbetter (mixed group with characteristics of both clusters), and none (not scoring high for genes enriched in either cluster) (Fig.7b). Both of the sequenced LGG tumors contained cells from all four groups in different proportions (Fig.7a, b). The signature top NMF genes of LGG were differentially expressed in cells in the four groups (Fig.7c). In a similar analysis of a pediatric single-cell dataset for NBL43 (Fig.7d), individual cells contained patterns of NBLpoor, NBLbetter, both NBLpoor and NBLbetter, and neither (Fig.7e, f), thus, supporting these epifactors as possible determinants of cellular states.

a UMAP plot of the LGG tumor cells color-coded (grey or black) based on the tumor sample from which they originate (LGG-03 or -04). b UMAP plot showing the assignment of tumor cells to four different groups based on the expression levels of the signature genes for the poor outcome (LGGlow) and better outcome (LGGhigh) tumor clusters for LGG. c Dot plot depicting the expression levels and the percent of cells expressing the signature genes for LGGlow and LGGhigh tumor clusters across the four groups determined using single-cell analysis. Only the genes that are differentially expressed across the four different groups are used for the plot. df Same plots as (ac), respectively, but for NBL tumor cells.

Read the original here:

Pan-cancer landscape of epigenetic factor expression predicts tumor outcome | Communications Biology - Nature.com