ARTICLE https://doi.org/10.1038/s41467-020-16000-6 OPEN MethCORR modelling of methylomes from formalin-fixed paraffin-embedded tissue enables characterization and prognostication of colorectal cancer Trine B. Mattesen 1, Mads H. Rasmussen 1, Juan Sandoval2,3, Halit Ongen 4, Sigrid S. Árnadóttir1, Josephine Gladov 1, Anna Martinez-Cardus 5,6, Manuel Castro de Moura7, Anders H. Madsen8, 1234567890():,; Søren Laurberg9, Emmanouil T. Dermitzakis 4, Manel Esteller 10,11,12,13, Claus L. Andersen 1,14 ✉ & Jesper B. Bramsen 1,14 ✉ Transcriptional characterization and classification has potential to resolve the inter-tumor heterogeneity of colorectal cancer and improve patient management. Yet, robust transcrip- tional profiling is difficult using formalin-fixed, paraffin-embedded (FFPE) samples, which complicates testing in clinical and archival material. We present MethCORR, an approach that allows uniform molecular characterization and classification of fresh-frozen and FFPE samples. MethCORR identifies genome-wide correlations between RNA expression and DNA methylation in fresh-frozen samples. This information is used to infer gene expression information in FFPE samples from their methylation profiles. MethCORR is here applied to methylation profiles from 877 fresh-frozen/FFPE samples and comparative analysis identifies the same two subtypes in four independent cohorts. Furthermore, subtype-specific prog- nostic biomarkers that better predicts relapse-free survival (HR = 2.66, 95%CI [1.67–4.22], P value < 0.001 (log-rank test)) than UICC tumor, node, metastasis (TNM) staging and microsatellite instability status are identified and validated using DNA methylation-specific PCR. The MethCORR approach is general, and may be similarly successful for other cancer types. 1 Department of Molecular Medicine, Aarhus University Hospital, 8200 Aarhus, Denmark. 2 Epigenomic Unit, Health Research Institute La Fe (ISSLaFe), Valencia, Spain. 3 Biomarker and precision medicine Unit, Health Research Institute La Fe (ISSLaFe), Valencia, Spain. 4 Genetic Medicine and Development, University of Geneva Medical School-CMU, 1 Rue Michel-Servet, 1211 Geneva, Switzerland. 5 Badalona Applied Research Group in Oncology (B-ARGO), Germans Trias i Pujol Research Institute (IGTP), Badalona, Barcelona, Catalonia, Spain. 6 Medical Oncology Service, Institute Catalan of Oncology (ICO), Badalona, Barcelona, Catalonia, Spain. 7 Josep Carreras Leukaemia Research Institute (IJC), Badalona, Barcelona, Catalonia, Spain. 8 Department of Surgery, Hospitalsenheden Vest, 7400 Herning, Denmark. 9 Colorectal Surgical Unit, Department of Surgery, Aarhus University Hospital, 8200 Aarhus, Denmark. 10 Josep Carreras Leukaemia Research Institute (IJC), Badalona, Barcelona, Catalonia, Spain. 11 Centro de Investigacion Biomedica en Red Cancer (CIBERONC), Madrid, Spain. 12 Institucio Catalana de Recerca i Estudis Avançats (ICREA), Barcelona, Catalonia, Spain. 13 Physiological Sciences Department, School of Medicine and Health Sciences, University of Barcelona (UB), Barcelona, Catalonia, Spain. 14These authors contributed equally: Claus L. Andersen, Jesper B. Bramsen ✉email: Cla@clin.au.dk; Bramsen@clin.au.dk NATURE COMMUNICATIONS | (2020)11:2025 | https://doi.org/10.1038/s41467-020-16000-6 | www.nature.com/naturecommunications 1 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-16000-6 C olorectal cancer (CRC) is a disease with extensive inter- analysis will identify far more associations and that these will patient heterogeneity, both molecularly and histopatho- include both functional gene-regulatory interactions and logically, which cannot be resolved by current clinical indirect associations e.g. between cell-type-specific RNA methods. Despite a continuous refinement of the UICC tumor, expression and cell-type-specific DNA methylation. We here node, metastasis (TNM) staging system to measure disease extent show that MethCORR, independent of whether the methylomes and define prognosis, the disease outcome still varies considerably were produced from FFPE or FF tissues, allows expression even for patients with the same tumor stage. Therefore, new information to be inferred for a large number of genes factors that can more precisely stratify patients into different risk (>11,000). Consequently, MethCORR enables a plethora of categories are clearly warranted1. molecular analyses to be performed on otherwise difficult-to- Recent attempts to resolve CRC heterogeneity and improve analyze FFPE tissues e.g. tumor characterization, tumor clas- prognostication include molecular subclassification and char- sification, and interpretation of expression signatures to derive acterization based on transcriptional profiling2–4. Consensus DNA methylation-based biomarkers. Hereby MethCORR also molecular subtype (CMS) classification stratifies CRC into four provides a path for improved, subtype-specific prognostication subtypes CMS 1–4 with distinct biology and histopathological of CRC using clinical FFPE samples. features2. Still, the CMS taxonomy itself has limited prognostic power for therapeutic decision-making5. To address this, we previously combined transcriptional subtyping with subtype- Results specific prognostic biomarkers to improve prognostication MethCORR infers RNA expression from DNA methylation. beyond TNM staging in retrospective cohorts3. This indicated a Here we developed the MethCORR approach that, by mapping clinical potential of using molecular classification and subtype- genome-wide correlations between RNA expression and DNA specific biomarkers as a complement to TNM staging for prog- methylation in FF samples, can infer gene expression information nostication. Furthermore, it highlighted the importance of in unrelated samples from their DNA methylation profiles. archived tumor material for biomarker discovery and pre-clinical Correlations were identified genome-wide using matching RNA validation. expression and 450K methylation data (methylation β-values) The strategies for transcriptional classification and subtype- from 394 FF CRC samples of The Cancer Genome Atlas (TCGA) specific prognostication were developed by, and still primarily Project, denoted the COREAD cohort (Fig. 1a and Supplementary rely on, profiling high-quality RNA purified from fresh-frozen Fig. 1a; Supplementary Table 1 and Supplementary Data 1). The (FF) tissue. However, high-quality RNA is often not recovered cohort was divided into two discovery sets (each n = 158) in from the formalin-fixed, paraffin-embedded (FFPE) tissue that is which genome-wide correlation analysis was performed inde- routinely archived in the clinic. This can preclude confident pendently and one validation set (n = 78; Fig. 1a). Our analysis transcriptional profiling and hereby complicate clinical testing of identified positively and negatively expression-correlated CpGs molecular classification and exploratory analysis in well-anno- (Spearman’s correlation P value < 0.01) overlapping in the two tated, archival FFPE material5–9. The clinical popularity of FFPE discovery sets for 17,776 of 20,530 genes (Fig. 1a). The majority of tissue is unlikely to change as it forms the basis for histopatho- the genes without expression-correlated CpGs were non- logical diagnoses and is a convenient, cost-effective preservation expressed (Supplementary Fig. 1b). To derive gene expression method. For wide utility, strategies for molecular classification, information for these 17,776 genes, we selected up to 200 CpGs characterization, and prognostication should therefore be com- whose methylation level were most negatively (≤100 sites) and patible with FFPE tissue. positively (≤100 sites) correlated with its expression (Fig. 1a). The Strategies based on DNA rather than RNA profiling may be a methylation levels of these expression-correlated CpGs were used way forward. DNA is considered less sensitive to degradation to calculate a MethCORR score (MCS) for each gene (formula in than RNA in FFPE samples10,11 and enzymatic strategies for Fig. 1b) and simple linear and polynomial regression modeling DNA repair have greatly improved the analysis of FFPE was used to identify genes with good correlations between MCSs DNA12–15. A strategy for robust analysis of clinical and archival and measured RNA expression (Fig. 1a). Models were established FFPE samples could involve DNA methylation as highly con- in the discovery sets by ten times tenfold cross validation and cordant DNA methylation profiles are produced from matched FF selected using root mean square error (RMSE) as a measure of and FFPE tissues when using the Illumina Infinium Human- model fit. We found good inter-sample correlations for 16,248 Methylation Beadchip technology14,16,17. In addition, many bio- genes in the discovery sets (R2 > 0.16) and confirmed these for logical traits, such as RNA expression and cell-type identity, are 11,222 genes in the validation set (gene model performances in associated with specific and robust DNA methylation patterns in Supplementary Data 2; Supplementary Fig. 1c–e). The 11,222 the genome18,19. This suggests that both gene expression and cell- genes were denoted MethCORR genes and the expression- type information may be extracted from DNA methylation pro- correlated CpGs of these define the COREAD MethCORR matrix files of FFPE samples and used for molecular classification and (≤200 CpGs × 11,222 genes; Supplementary Data 3) that was used prognostication, as an alternative to RNA profiling. Furthermore, for calculation of MCSs from DNA methylation profiles of all conversion of methylation profiles into a gene-centric expression samples analyzed in this study (Fig. 1c). We also investigated if format would allow molecular analysis of FF and FFPE samples RNA expression was better modeled using the ≤200 expression- using the plethora of bioinformatics tools, databases, and sig- associated CpGs for each gene directly, instead of using MCSs, natures established for RNA expression analysis. but found no improvement in overall performance (R2 and Motivated by this, we here develop MethCORR, an approach, RMSE; Supplementary Fig. 1f). Similarly, adding age and gender which identifies genome-wide correlations between gene information to MCS-based models did not improve overall per- expression and DNA methylation and use this to obtain gene formances (Supplementary Fig. 1g). This likely reflect that CRC- expression and cell-type information in independent samples induced methylation changes are much greater than the subtler from their DNA methylation profiles. In homogenous cell effects of age and gender in normal tissues23. Still, MethCORR preparations, associations between gene expression and DNA captures gender-specific expression by including CpGs located on methylation have been observed only for a small fraction of chromosome X and Y in the MethCORR matrix. Accordingly, genes when analyzing local promoters, gene bodies, or nearby known gender-specific RNAs exhibited gender-specific inferred enhancers20–22. We hypothesize that genome-wide correlation RNA expression (Supplementary Fig. 1h). 2 NATURE COMMUNICATIONS | (2020)11:2025 | https://doi.org/10.1038/s41467-020-16000-6 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-16000-6 ARTICLE a TCGA COREAD c Genes expression correlated CpG sites MethCORR Data analysis correlation matrix Set 1 Set 2 Set 3 n = 158 n = 158 n = 78 FF or FFPE sample Calculate MCSs Use MCSs to perform RNA expression- with DNA methylation for 11,222 RNAs using a analysis of interest DNA methylation data (450K/EPIC) pre-established MethCORR correlation for matrix 20,530 RNAs Common top ranking CpG sites for 17,776 RNAs d 20 R 2 = 0.93 e COREAD SYSCOL f 100% matrix matrix FF tissue (COREAD, pt. 6882) Fraction of samples RMSE=0.95 1 Intra-sample R 2,RNA vs iRNA 80% Calculation of MCSs for each iRNA expression (log2) MCSs for 15 RNA from top ranking CpG sites 17,776 RNAs 0.8 60% cross-validation 10 0.6 40% 10 x 10-fold Model relationship between 20% each RNA-specific MCS and 5 0.4 RNA expression 0% 0.2 CMS1 CMS2 CMS3 CMS4 0 RNA expr. RNA expr. RNA expr. RNA expr. Select model with 0 5 10 15 20 (n = 61) (n = 198) (n = 40) (n = 71) 0 lowest RMSE RNA expression (log2) C CMS1 iRNA CMS3 iRNA SY C = SY L (n O O 78 (n FF tissue (COREAD, pt. 6882) SC set R SC R ) EA t 3 EA CMS2 iRNA CMS4 iRNA O O D D L se (n 17,275 simple linear = = Independent 3 31 39 (n regression models 4) 4) = validation of 62 501 polynomial ) selected models regression models Exclude models with MCS-RNA correlation ≤ 0.4 g 12 R 2 = 0.69 FFPE tissue (COREAD, pt. 6781) 12 FFPE tissue (COREAD, pt. 6781) (R 2 ≤ 0.16) in discovery and validation sets R 2 = 0.91 RMSE = 1.00 10 RMSE = 0.57 10 iRNA expression (log2) RNA expression (log2) 8 8 11,222 MethCORR genes with well-performing 6 6 regresssion models 4 4 2 2 0 0 b 0 2 4 6 8 10 12 0 2 4 6 8 10 12 RNA expression (log2) RNA expression (log2) FF tissue (COREAD, pt. 6781) FF tissue (COREAD, pt. 6781) Fig. 1 Development of the MethCORR approach, MethCORR scores, and inferring of RNA expression. a Overview of MethCORR development in the COREAD cohort using matched RNA-sequencing and 450K methylation data. The cohort was divided into two discovery sets (each n = 158) and one validation set (n = 78). Genome-wide RNA expression-DNA methylation correlations were identified in each discovery set and shared top expression- correlated CpGs for each RNA were selected (≤100 positively and ≤100 negatively correlated CpGs; Spearman’s correlation P value < 0.01). A MCS was calculated for each gene using DNA methylation β-values of expression-correlated CpGs and the formula given in (b). RNA expression of each gene was modeled from its MCS using simple linear- and polynomial regression models and 10 × 10-fold cross validation in set 1 + 2. Simple linear models were selected for all RNAs except when polynomial models exhibited a ≥5% decrease in RMSE values (Supplementary Fig. 1d). Only models with R2 > 0.16 between inferred RNA (iRNA) and observed RNA expression in both the discovery set and the independent validation set 3 was kept for further analysis (n = 11,222, termed MethCORR genes). b Formula for calculating MCSs from DNA methylation β-values. c Overview of MethCORR applications. Fresh- frozen (FF)/FFPE CRC samples with 450K/EPIC methylation profiles can be applied to the MethCORR matrix for calculation of MCSs and iRNA expression. d Scatterplot showing intra-sample correlation between iRNA and RNA expression in a representative COREAD validation sample. e Plot showing R2 of intra-sample iRNA and RNA expression correlations for all samples of the COREAD validation set 3 and SYSCOL cohort when using the COREAD-derived MethCORR matrix (left) and for all samples of the SYSCOL validation set 3 and COREAD cohort when using the SYSCOL-derived MethCORR matrix (right). f Histogram showing overlap in CMS subtype predictions in COREAD CRC samples using RNA expression or iRNA expression for classification. g Scatterplot showing intra-sample correlation between iRNA (left) or RNA expression (right) from a FFPE sample and RNA expression in a matched fresh- frozen COREAD sample. Next, we investigated characteristics of the MethCORR genes Data 4) and SYSCOL samples (median R2 = 0.76 (range = included in the MethCORR matrix. MethCORR genes exhibited 0.62–0.82); Fig. 1d–e; Supplementary Data 5). To evaluate the greater variation in RNA expression (Supplementary Fig. 2a), were robustness of MethCORR to differences between cohorts, we more frequently dysregulated in cancer vs. normal mucosa repeated the entire MethCORR discovery and validation process (Supplementary Fig. 2b) and encompassed relatively fewer house- using the SYSCOL cohort to construct a SYSCOL MethCORR hold genes (Supplementary Fig. 2c) than the set of genes not matrix, derive MCSs, and to infer iRNA expression (Fig. 1a; included in the MethCORR matrix. Importantly, the MethCORR Supplementary Data 6–7). Again, we found high intra-sample genes exhibited the same stroma score distribution as the full set of correlations between observed RNA and iRNA expression (SYSCOL genes (Supplementary Fig. 2d). This indicates that MethCORR set 3, median R2 = 0.92 (range = 0.87–0.95); COREAD median maintains the ability to characterize both traits of the cancer cells R2 = 0.74 (range = 0.55–0.82); Fig. 1e; Supplementary Data 4 and and the surrounding stroma. The established MCS regression 5). We speculated that the moderate decrease in R2 between cohorts models were next used to calculate inferred RNA (iRNA) expression was caused by differences in RNA quantification methods rather for MethCORR genes in the validation samples of the COREAD than the MethCORR approach. In support, comparative analysis cohort (set 3) and in an independent Danish CRC cohort, denoted of COREAD validation samples using normalized RNA expres- SYSCOL3. We found a high intra-sample correlation between sion data from the UCSC XENA database24 and the National measured RNA and iRNA expression in the COREAD validation Cancer Institute (NCI) genomic database commons (GDC)25 samples (median R2 = 0.93 (range = 0.82–0.96); Supplementary confirmed that MethCORR iRNA-RNA correlations were not NATURE COMMUNICATIONS | (2020)11:2025 | https://doi.org/10.1038/s41467-020-16000-6 | www.nature.com/naturecommunications 3 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-16000-6 Table 1 R2 and RMSE for intra-sample correlations between MethCORR inferred RNA expression (iRNA), RNA expression, or MCS in FFPE samples and RNA expression or MCS in matched fresh-frozen tissue. TCGA COREAD R2 iRNA (FFPE) R2 RNA (FFPE) vs. R2 MCS (FFPE) vs. RMSE iRNA (FFPE) RMSE RNA (FFPE) RMSE MCS (FFPE) patient Id vs. RNA (FF) RNA (FF) MCS (FF) vs. RNA (FF) vs. RNA (FF) vs. MCS (FF) Pt. 6650 0.94 0.87 1.00 0.47 0.69 0.04 Pt. 5659 0.92 0.74 1.00 0.54 1.08 0.03 Pt. 5661 0.92 0.67 0.99 0.54 1.25 0.03 Pt. 5665 0.91 0.72 0.98 0.57 1.02 0.04 Pt. 6781 0.91 0.69 0.98 0.54 1.00 0.03 Pt. 6780 0.90 0.81 0.99 0.60 0.82 0.03 Pt. 2684 0.88 0.67 0.98 0.65 1.03 0.04 Pt. 3810 0.87 0.70 1.00 0.66 0.98 0.02 Pt. 5656 0.80 0.63 0.98 0.83 1.11 0.07 lower than if applying two different RNA normalization strategies consensus clustering was performed in TNM stage II–III to the same samples (Supplementary Fig. 2e). COREAD and SYSCOL samples with available 450K methylation In accordance with the high intra-sample correlations between data and in two independent FFPE TNM stage II–III cohorts, measured RNA and iRNA expression, we found a good overlap in denoted FFPE1 and FFPE2 (Supplementary Table 1 and Sup- CMS (84% agreement) and CRC intrinsic subtype (CRIS; 75% plementary Data 12). Our focus was on stage II–III patients, agreement) predictions when using the measured RNA or iRNA which are most relevant for prognostic biomarker identification expression as input (Fig. 1f and Supplementary Fig. 2f). due to their heterogeneous prognosis1. Two MethCORR subtypes, In situations where high-quality RNA is not obtainable, iRNA CRC1 and CRC2, were identified in all four cohorts (Supple- expression may provide better estimates of gene expression than mentary Fig. 3c) and Submap analysis28 confirmed the corre- RNA sequencing, as even moderate declines in RNA quality can spondence between the CRC1 and CRC2 subtypes in the different lead to unreliable expression profiles26,27. Indeed, samples with cohorts (Supplementary Fig. 3d; FDR < 0.05). In agreement, the lowest correlation between measured RNA and iRNA samples clustered according to subtype in a PCA of all four CRC expression had significantly lower RNA quality than high cohorts together, irrespectively of their preservation-type status correlation samples (P value < 0.0001, Wilcoxon rank sum (Supplementary Fig. 3e). We next performed comparative sub- (WRS) test; Supplementary Fig. 2g). In contrast, no equivalent type characterization in all cohorts, which indicated that CRC1 drop in 450K methylation data quality was observed (Supple- and CRC2 differed in terms of DNA methylation, chromosomal mentary Fig. 2g). Compromised RNA quality is inherent to FFPE instability, and stromal/immune cell activity (Fig. 2a and Sup- tissue10,11. In agreement, analysis of nine COREAD samples with plementary Fig. 3f). These are well-known characteristics for the available RNA sequencing and 450K methylation profiles from serrated/microsatellite instability status (MSI) and conventional matched FF and FFPE tissues identified higher intra-sample R2’s CRC pathways, respectively, pointing to a biological relevance of between FF RNA sequencing and FFPE iRNA profiles (median the MethCORR subtypes. R2 = 0.91 (range: 0.80–0.94)) than between FF and FFPE RNA- Further subtype characterization was performed using pre- sequencing profiles (median R2 = 0.7 (range: 0.63–0.87); P ranked gene set enrichment analysis (GSEA)29. Initially, we value < 0.001, WRS test; Fig. 1g and Table 1; Supplementary investigated if similar gene set enrichments were identified when Data 8–11 and Supplementary Table 2). MCSs from matched using MCSs vs. RNA expression as input (Fig. 2b) or when MCSs FFPE and FF samples were even higher correlated (median R2 = were derived from FF vs. FFPE samples (Fig. 2c). Indeed, a high 0.98 (range: 0.98–1.00); Table 1), which likely reflect that 450K concordance was observed between normalized enrichment methylation profiles were themselves highly correlated (median scores for most gene sets in both situations, supporting that R2 = 0.96 (range: 0.94–0.98); Supplementary Fig. 2h), as reported expression-correlated MCSs can substitute RNA expression and previously14,16,17. Additional evidence came from principal enable analysis of FFPE tissue. MCS-based GSEA of each cohort component analysis (PCA). Here samples clustered according uniformly showed that the CRC1 subtype was enriched in gene to preservation method when analyzing FF and FFPE RNA- sets associated with immune- and stromal processes/cell types sequencing profiles together, whereas samples clustered more such as inflammation, epithelial-mesenchymal transition (EMT), according to patient ID when analyzing RNA profiles of FF cancer-associated fibroblasts (CAFs), and T/B cells (Fig. 2d and samples together with iRNA or MCS profiles of FFPE samples Supplementary Table 3). Furthermore, CRC1 was enriched in (Supplementary Fig. 2i). gene sets associated with positive MSI-, CIMP-, and serrated Collectively, this showed that MethCORR expression measures CRC-status, whereas CRC2 tumors were enriched in gene sets (MCSs and iRNAs) can be inferred from DNA methylation for a associated with conventional CRC and a more undifferentiated large number of genes, even when methylation data are based on cell status (Fig. 2d and Supplementary Table 3). Similar results FFPE tissue. were obtained for the two FF cohorts when using RNA expression as input, rather than MCSs (Fig. 2d). Despite biological MethCORR identifies two subtypes in FF and FFPE cohorts. differences, no difference in relapse-free survival (RFS) was We next investigated if inferred expression profiles allow uniform observed between CRC1 and CRC2 (Fig. 2e). subtype discovery and characterization of both FF and FFPE Collectively, these results demonstrate that MethCORR allows cohorts using bioinformatics strategies normally reserved for FF uniform discovery and characterization of biologically relevant samples with high-quality RNA expression profiles. As input, we CRC subtypes in FF and FFPE samples using well-established employed MCS profiles as they strengthen the focus on cancer bioinformatics tools. cell-related traits during subtype discovery as compared with RNA and iRNA profiles (Supplementary Fig. 3a, b). Subtype A MethCORR map characterizes CRC subtypes. By analysis of discovery by non-negative matrix factorization (NMF)-based expression-correlated CpGs in the MethCORR matrix, we found 4 NATURE COMMUNICATIONS | (2020)11:2025 | https://doi.org/10.1038/s41467-020-16000-6 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-16000-6 ARTICLE a COREAD FF1 FFPE1 CRC1 (n = 48) CRC2 (n = 86) CRC1 (n = 51) CRC2 (n = 53) MethCORR CRC subtypes MSI MSS NA CIN Score Low High DNA methylation score Low High ESTIMATE Stroma Score Low High ESTIMATE Immune Score Low High b CRC1 vs. CRC2 d CRC1 vs CRC2 CRC1 vs CRC2 with MethCORR scores with RNA expr. NES using RNA expression 6 Pearson’s r = 0.91 GSEA gene set FF1 FF2 FFPE1 FFPE2 FF1 FF2 4 CAFs 3.27 2.49 2.54 2.96 2.53 0.76 P < 10 –100 in COREAD FF1 Endothelium 3.19 2.41 2.65 3.09 2.62 1.41 2 EMT 2.96 2.35 2.46 2.55 2.66 2.09 Tumor microenvironment 0 T-cells 3.96 3.45 3.14 3.00 2.76 2.90 -8 –6 –4 –2 0 2 4 6 8 B-cells 3.02 2.40 2.34 2.33 2.39 2.19 –2 Macrophages 2.78 2.12 2.21 2.38 2.04 1.18 Significant both (n = 4741) Dendritic cells 3.15 2.59 2.63 3.00 2.52 1.95 –4 Significant in one (n = 2163) Myeloid derived suppressor cells 2.89 2.34 2.23 2.57 2.04 1.15 –6 Hallmark inflammatory response 3.72 3.10 3.19 3.30 2.61 2.28 Non-significant both (n = 7106) Hallmark IFNg response 3.85 3.49 3.41 3.35 2.58 2.95 –8 Hallmark complement 2.78 2.44 2.44 2.58 2.14 2.01 Enriched NES using MCSs in COREAD FF1 Hallmark allograft rejection 3.67 3.06 2.83 2.83 2.76 2.46 in CRC1 Wound healing 2.32 1.78 2.02 2.36 2.04 1.48 4.4 NES c CRC1 vs. CRC2 Up in CIMP vs. non-CIMP CRC Up in MSI vs. MSSCRC 2.99 1.59 4.43 2.76 4.43 2.85 2.76 1.30 2.36 1.68 3.75 3.00 0 –4.9 NES using MCSs in FFPE1 6 Up in BRAF mut. vs. non-BRAF mut. CRC 2.50 3.59 3.74 2.20 2.09 3.25 Pearson’s r = 0.90 Up in serrated vs. conventional CRC 2.54 2.90 2.86 2.34 2.20 2.92 Enriched Cancer epithelial cell processes P < 10–100 4 in CRC2 Genes upreg. by APC 2.40 2.11 2.05 2.05 1.98 1.50 2 8 Up in colon crypt top vs. bottom 1.56 1.73 1.64 1.74 1.38 1.24 Up in non-CIMP vs. CIMP CRC –4.41 –4.59 –4.61 –3.73 –4.15 –3.67 0 Up in MSS vs. MSI CRC –3.75 –3.41 –3.61 –3.42 –3.36 –2.82 –8 –6 –4 –2 0 2 4 6 –2 Up in non-BRAF mut. vs. BRAF mut. CRC –4.87 –4.14 –4.35 –4.12 –4.29 –3.15 Significant both (n = 3912) Up in conventional vs. serrated CRC –3.84 –3.45 –3.36 –3.63 –3.37 –2.69 –4 Genes downreg. by APC –2.29 –1.57 –1.78 –2.28 –1.82 –1.27 Significant in one (n = 2546) Up incolon crypt bottom vs. Top –2.56 –1.96 –2.14 –2.53 –1.52 –1.43 –6 Non-significant both (n = 7552) Up in WNT pathway –2.04 –1.78 –1.80 –1.91 –1.38 –1.45 –8 Undifferentiated cancer –2.78 –1.91 –2.26 –2.74 –1.78 –1.29 NES using MCSs in COREAD FF1 Hallmark MYC targets –4.43 –3.20 –3.52 –4.52 –3.19 –1.10 Hallmark DNA repair –2.28 –1.77 –1.86 –2.51 –1.87 –1.43 Hallmark E2F targets –3.93 –2.44 –2.90 –3.96 –2.68 –0.80 Hallmark G2M checkpoint –3.25 –1.98 –2.26 –3.28 –2.04 –0.84 e FF1-FF2 patients stratified by FFPE1 patients stratified by MethCORR subtypes MethCORR subtypes 1.00 1.00 Proportion patients Proportion patients w/o recurrence w/o recurrence 0.75 0.75 0.50 0.50 0.25 0.25 0.00 NS 0.00 NS 0 50 100 150 0 10 20 30 40 Time (Months) Time (Months) Number at risk Number at risk CRC1 123 66 20 1 CRC1 51 46 38 37 0 CRC2 164 91 28 0 CRC2 53 50 41 38 0 CRC1 CRC2 CRC1 CRC2 Fig. 2 MethCORR based NMF clustering identifies the same two CRC subtypes in fresh-frozen and FFPE cohorts. a Main molecular features of the CRC1 and CRC2 MethCORR subtypes in the COREAD FF1 and the FFPE1 cohort (Supplementary Table 1). MSI and MSS status is indicated in black and white. CIN scores were derived for COREAD and FFPE1 samples using GISTIC and EPIC DNA methylation data, respectively, and sample DNA methylation scores were calculated as the 40th percentile of DNA methylation β-values for all CpGs. Stroma- and Immune Scores were generated from MCSs using the ESTIMATE software69. b–c Scatterplots showing the correlation between normalized enrichment scores (NESs) for ~17 K gene sets of The Molecular Signatures Database (MSigDB) v6.1 from a pre-ranked GSEA of CRC1 vs. CRC2 subtypes in the COREAD FF1 cohort using either MCSs (X-axis) or RNA expression (Y-axis) as input (b) and a pre-ranked GSEA of CRC1 vs. CRC2 subtypes in either the COREAD FF1 cohort (X-axis) or FFPE1 cohort (Y-axis) using MCSs as input (c). Pearson’s r and P value (Wilcoxon rank sum test) is indicated. d Table showing selected gene sets differentially enriched between CRC1 and CRC2 subtypes as evaluated by pre-ranked GSEA performed using MCSs or RNA expression in the fresh-frozen COREAD FF1 and SYSCOL FF2 cohorts and MCSs for the FFPE cohorts (Supplementary Table 1). Gene sets with positive NES are enriched in CRC1 (red colors), whereas negative NES indicate enrichment in CRC2 (blue colors). Gene sets enriched/depleted at a high significance level are highlighted in bold (FDR < 0.05). See methods section and Supplementary Table 3 for origin of gene sets. e Kaplan–Meier plots showing the relapse-free survival of CRC patients stratified according to subtype. Left panel: patients with fresh-frozen tumors and good clinical follow-up (the COREAD FF1 and SYSCOL FF2 cohorts; Supplementary Table 1) were combined to increase the number of relapse events. Right panel: patients with FFPE tumors and good clinical follow-up (The FFPE1 cohort; Supplementary Table 1). Significance was evaluated by the log-rank test. NATURE COMMUNICATIONS | (2020)11:2025 | https://doi.org/10.1038/s41467-020-16000-6 | www.nature.com/naturecommunications 5 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-16000-6 that most CpGs were not located on the same chromosome as the overlain with MethCORR interpretations of RNA-based bio- gene they correlate with (Supplementary Fig. 4a). Instead, the markers and signatures defining CAFs, endothelium, myeloid most frequently occurring CpGs were located in genomic regions cells, T cells, and B cells (Supplementary Fig. 5e). Hence, the that exhibited great cell-type-specific variation in DNA methy- MethCORR map can suggest cell types associated with RNA lation, as evaluated in 17 tissue types (GSE5019218; Supplemen- biomarkers and signatures via comparison to known cell-type- tary Fig. 4b). Hence, the MethCORR matrix may help associate specific methylation profiles. gene expression with particular cell types by comparing the Based on this, we envisioned that the MethCORR map would methylation pattern of expression-correlated CpGs to known visualize and suggest sources of inter-tumor heterogeneity DNA methylation (or DNAse I hypersensitivity) profiles of cell between and within subtypes. CRC heterogeneity can arise from monocultures/homogenous cell preparations. Indeed, expression- both differences in TME cell composition and in the differentia- correlated CpGs for the T-cell-specific CD3 Epsilon (CD3E) gene tion status of tumor epithelial cells. For example compared with overlapped with T-cell specific DNase I hypersensitive sites and normal mucosa, CRCs can lose mature enterocyte traits and DNA methylation patterns characteristic of T-cells (Supplemen- rather resemble enterocyte precursors, transit amplifying (TA) tary Fig. 4c, d). Similarly, expression-correlated CpGs for fibro- and stem cells, or undergo EMT2,33,34. Mapping of MCS z-scores blast activation protein alpha (FAP) and epithelial cellular from individual tumors revealed inter-tumor heterogeneity in adhesion molecule (EPCAM) overlapped with patterns char- both subtypes. For CRC1, heterogeneity was pronounced in the acteristic of stromal cells/fibroblasts and intestinal epithelial cells, TME cluster and few samples had a dominant epithelial pattern respectively (Supplementary Fig. 4c, d). We also found that the (Fig. 3f). Three TME patterns were frequently observed, one genes with greatest expression-correlated CpG site overlap with overlapping with CAF/fibroblast (CAF/fibroblast pattern), CD3E, FAP, and EPCAM were themselves significantly associated another with CD14+ monocytic cells/platelets (inflammation with T-, stromal/fibroblast-, and epithelial-cell activities as eval- pattern), and the last with lymphocytic T cells and B cells uated by gene list enrichment analysis30 (Supplementary Fig. 4e; (lymphocyte pattern; Fig. 3e–g). This suggested that TME cell P value < 0.05 by the Enrichr software30). This showed that composition is a major contributor to intra-subtype heterogeneity analysis of expression-correlated CpGs help identify clusters of in the immune-infiltrated CRC1 subtype. The TME patterns were co-expressed genes and link them to particular cell types via less dominant among CRC2 samples (Fig. 3h) consistent with comparison to cell-type-specific DNA methylation profiles. CRC2 conventional-like tumors being less immune-infiltrated2 To analyze expression correlations in a genome-wide format, (Fig. 2a, d). Instead, CRC2 heterogeneity was pronounced within we created a MethCORR map by clustering all MethCORR genes epithelial map regions and four patterns were observed (Fig. 3h): according to their overlap in expression-correlated CpGs (Fig. 3a). Two regions were dominated by signatures of enterocyte Foremost, the map was used to visualize differences between precursors and TA cells as estimated by overlapping the map CRC1 and CRC2 by coloring gene nodes according to their with RNA signatures defining specific differentiation states of difference in median MCS z-score between the subtypes intestinal epithelial cells33 (Fig. 3i). A third region overlapped (Δmedian z-score; Fig. 3a). The differences were near-identical with a mature enterocyte signature characteristic of normal for FF and FFPE cohorts (Fig. 3a, b and Supplementary Fig. 5a; mucosa samples (Fig. 3i and Supplementary Fig. 5f). Finally, an Δmedian z-score Pearson’s r range: 0.88–0.97, P value < 10−100, EMT pattern was identified in CRC2 by overlaying the map with WRS test) and near-identical to a MethCORR map comparing MCSs of Hela cells undergoing EMT35 (Fig. 3i) and GSEA serrated/MSI and conventional adenocarcinomas from the 450K showed enrichment of EMT signatures in the CRC2 samples with methylation dataset GSE6806031 (Fig. 3c; Δmedian z-score this EMT pattern (as compared with an early enterocyte pattern; Pearson’s r range: 0.87–94, P value < 10−100, WRS test). Similar Supplementary Fig. 5g). Collectively, this suggested that epithelial results were obtained when the map was overlain with differentiation status is an important contributor to heterogeneity MethCORR interpretation of a transcriptional gene set defining in the CRC2 subtype. Finally, the above heterogeneity was also serrated vs. conventional CRC (Supplementary Fig. 5b; Pearson’s identifiable among CRC cell lines and CMS subtypes (Supple- r range = 0.94–98, P value < 10−100, WRS test; for comparison to mentary Fig. 5h, i). MSI status, CIMP status, CMS- and CRIS-classification status see Supplementary Fig. 5c, d). This suggested that CRC1 and CRC2 subtypes resembles serrated/MSI and conventional carci- MethCORR interprets prognostic RNA signatures. We next nomas, respectively. In support, Submap analysis confirmed that investigated if MethCORR would also help identify DNA CRC1 and CRC2 subtypes from all four cohorts corresponded to methylation-based biomarkers suited for prognostication using the serrated/MSI and conventional subtypes from the GSE68060 FF and FFPE samples. Our strategy was to use the MethCORR dataset31 (Supplementary Fig. 3d). Furthermore, CRC2 encom- map to interpret established, prognostic RNA signatures and passed several map regions associated with high CIN scores, suggest cell types associated with tumor aggressiveness, which can whereas CRC1 encompassed a large tumor microenvironment be evaluated in DNA samples based on the cell-type specificity of (TME) cluster characterized by genes with high stroma scores, as methylation. Analysis of five prognostic signatures, CRC-11336, expected for conventional and serrated/MSI tumor subtypes2,32, ColoGuideEx37, Oncotype DX38, ColoPrint39, and Tian et al.40 respectively (Fig. 3d). showed that MCSs for almost all stromal transcripts were posi- tively correlated with the median MCS for all signatures (Fig. 4a). This suggested that all signatures associated high TME activity The MethCORR map characterizes intra-tumor heterogeneity. with poor prognosis. MethCORR map analysis of the signatures To investigate the large TME cluster in greater detail and provide revealed two distinct patterns within the TME cluster: The CRC- insight into sources of CRC heterogeneity, the map was overlain 113, ColoGuideEX, and Oncotype DX signatures associated with with MCS z-scores calculated from DNA methylation profiles of a CAF-like pattern (Figs. 3e, f, and 4b), cancer invasiveness and epithelial, immune, stem, and mesenchymal cells (primarily cell hepatocyte growth factor (HGF) expression41 (Fig. 4c, d). The monocultures; Supplementary Table 4 and Supplementary ColoPrint and Tian et al. signatures (Fig. 4e) associated with Data 13). This identified map regions representing CAFs, CD14+ an inflammation/wound healing pattern (Figs. 3e, f, and 4c) monocytes, CD3+ T cells, and CD19+ B cells among others encompassing blood platelets, CD14+ monocytes (Fig. 3e), and (Fig. 3e). Again, similar results were obtained when the map was transforming growth factor beta 1 (TGFB1) expression (Fig. 4d). 6 NATURE COMMUNICATIONS | (2020)11:2025 | https://doi.org/10.1038/s41467-020-16000-6 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-16000-6 ARTICLE a b c CRC1 Serrated/MSI enriched CRC enriched MCS z-score MCS z-score ΔMedian ΔMedian 1.0 1.0 –1.0 –1.0 CRC2 Conventional enriched CRC enriched CRC1 d e enriched 20p11–13 MCS z-score ΔMedian 1.0 8p21–22 –1.0 CRC2 enriched 8q21–24 20q KRT7 (Serrated CRC marker) ALPI (mature enterocyte marker) 11–13 13q12–14 ANXA10 (Serrated CRC marker) CHGA (Enteroendicrine cell marker) TME cluster KRT23 (MSS CRC marker) FCGBP (Goblet cell marker) 7p13–22 CD3+ T-cells CD19+ B-cells CTNNB1 (Conv. CRC marker) CD14+ Monocytes Platelets (RNA profile) High Stroma score High CIN score Endothelium Stem cells (embryonic) CAFs f g CRC1: Epithelial pattern CRC1: CAF/Fibroblast pattern CRC1: Inflammation pattern CRC1: Lymphocyte pattern Cytotoxic T-cells Leukocytes Hallmark inflammatory MCS z-score MCS z-score MCS z-score MCS z-score 2 2 2 2 response CAFs –2 Pt. 5661 –2 Pt. 6463 –2 Pt. 6304 –2 Pt. 6895 h i CRC2: Early enterocyte precur./TA CRC2: Late precur./mat. enterocyte CRC2: Mature enterocyte CRC2: Epithelial EMT TA cells Early enterocyte precursor Late enterocyte precusor MCS z-score MCS z-score MCS z-score MCS z-score 2 2 2 2 Enterocytes Epithelial cell –2 Pt. a1dd –2 Pt. a28c –2 Pt. 6136 –2 Pt. 6460 EMT Fig. 3 A MethCORR map identifies characteristics of CRC subtypes and intra-subtype heterogeneity. a The MethCORR map is a representation of the MethCORR matrix established by clustering genes (cluster nodes) according to their overlap in expression-correlated CpGs (cluster edges) using Enrichment Map63. Each gene is colored according to the difference in median MCS z-scores (ΔMedian MCS z-score) comparing CRC1 and CRC2 within the COREAD FF1 cohort (Supplementary Table 1). Epithelial and CRC-related genes are highlighted by circles. b MethCORR map with genes colored according to ΔMedian MCS z-scores comparing CRC1 and CRC2 within the FFPE1 cohort. c MethCORR map with genes colored according to ΔMedian MCS z-scores comparing serrated/MSI and conventional CRCs (GSE6806031). d MethCORR map with genes colored according to a high stroma score (≥0.548; red) or high CIN score (≥0.4; blue). A cluster encompassing genes with high stroma scores was named tumor microenvironment (TME) cluster. e Magnification of the TME cluster with genes colored according to high MCS z-scores for either CD3+ T cells, CD19+ B cells, CD14+ monocytes, platelets (RNA profile; MSigDB M7732), endothelium, stem cells (embryonic), or CAFs. MCS z-score profiles were calculated within a set of public DNA methylation profiles of cell monocultures and tissues (Supplementary Table 4 and Supplementary Data 13). Black lines indicate separation of the TME into lymphocyte, inflammation, and CAF/stem cell regions based on differences in cell-type composition. f and h MethCORR maps with genes colored according to the MCS z-scores of representative CRC1 (f) and CRC2 (h) samples calculated within all samples of the COREAD FF1 cohort. Black lines indicate TME patterns. g MethCORR map with genes colored according to high correlation to median MCS (cMCS) for three gene sets defining either cytotoxic T cells (MSigDB M1324729 (BioCarta)), leukocytes48, hallmark inflammatory response29, or CAFs48). i MethCORR map with genes colored according to high cMCS for transcriptional gene sets up in either transit amplifying (TA) cells, early enterocyte precursors, late enterocyte precursors, or enterocytes33. Genes with >5% increase in MCSs during EMT of epithelial HeLa cells35 are indicated in red. See methods section for details of ΔMedian MCS z-score and cMCS calculations. Hence, the prognostic signatures overlapped in predictions, and signatures denoted SSC prognosis and CIN prognosis. These are pointed to CAF or inflammation/wound healing as associated intended for immune-infiltrated/serrated and conventional car- with poor prognosis CRC. We recently reported that subtype- cinoma subtypes3, which correspond to CRC1 and CRC2 in this specific RNA signatures can improve prognostication beyond study, respectively. MethCORR map analysis suggested that TNM staging in multiple CRC cohorts3. Therefore, MethCORR depletion of immune cells, including T cells, was associated with was also used to interpret these subtype-specific prognostic the SSC prognosis signature (Figs. 3e and 4c, f), whereas a CAF NATURE COMMUNICATIONS | (2020)11:2025 | https://doi.org/10.1038/s41467-020-16000-6 | www.nature.com/naturecommunications 7 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-16000-6 a b 100% Fraction of stromal median MCS (cMCS) median MCS (cMCS) median MCS (cMCS) Correl. to signature Correl. to signature Correl. to signature 75% transcripts 1 0.3 1 50% 25% 0% –1 –0.3 –1 CCR113 ColoGuideEX Oncotype DX Pos. correl. to signature median MCS Neg. correl. to signature median MCS c d e median MCS (cMCS) median MCS (cMCS) Correl. to signature Correl. to signature 0.75 1 –0.75 –1 Coloprint Tian et al. High cMCS to cytotoxic T-cells High correl. to IL2 MCS High cMCS to wound healing High correl. to TGFB1 MCS High cMCS to cancer invasiveness High correl. to HGF MCS f g h MCS % increase MCS % increase Correl. to signature Correl. to signature median ΔcMCS median ΔcMCS 0.25 0.25 0.15 0.30 –0.25 –0.25 –0.15 –0.30 COREAD SYSCOL SSC prognosis signature CIN prognosis signature TNM st. IV vs I TNM st. IV vs I i j FF1-FF2 FFPE1 0.1 0.3 0.08 Correlation to recurrence status Correlation to recurrence status 0.2 0.06 0.1 0.04 0.02 0 0 –0.1 –0.02 –0.2 –0.04 –0.3 –0.06 High correl. to ACTA2 MCS High correl. to CD3E MCS –0.08 –0.4 CRC1 CRC2 CRC1 CRC2 High correl. to PDPN MCS High correl. to HNF4A MCS ACTA2/PDPN (Fibroblast/EMT) top50 correlated genes CD3E (T-cell) top50 correlated genes Fig. 4 The MethCORR map suggests cell types associated with prognostic RNA signatures. a Bar plot showing the fraction of stromal genes (stroma score > 0.5) that have positive or negative cMCSs (correlation to the median MCS) calculated for five prognostic RNA signatures CRC-11336, CologuideEx37, Oncotype DX38, Coloprint39, and Tian et al.40 (see “Methods” for calculation of cMCSs). Stromal transcripts were significantly enriched among positively vs. negatively correlated transcripts for all five signatures (P value < 10−100, Wilcoxon rank sum test). b Magnification of the TME cluster, where genes with the highest cMCSs for the prognostic CRC-11336 (left), CologuideEx37 (middle), and Oncotype DX38 (right) signatures are highlighted. c Magnification of the TME cluster, where genes with the highest cMCSs for published gene sets defining cytotoxic T-cells (MSigDB M1324729 (BioCarta); green), wound healing (MSigDB M1195729,73; red), or cancer invasiveness (MSigDB M257229,74; orange) are highlighted. d Magnification of the TME cluster, where genes with the highest correlation to the MCS of the IL2 (green), TGFB1 (red), and HGF (orange) genes are highlighted. e Magnification of the TME cluster, where genes with the highest cMCSs for the prognostic Coloprint39 (right) and Tian et al.40 (left) signatures are highlighted. f Magnification of the TME cluster, where genes with the highest ΔcMCSs for the prognostic SSC prognosis signature3 and g the CIN prognosis signatures3 are highlighted. h Magnification of the TME cluster colored according to the gene-specific percentage change in median MCSs between TNM stage I and IV CRCs of the COREAD cohort (left) and SYSCOL cohort (right). i Magnification of the TME cluster, where genes with the highest correlation to the MCS of CD3E (green), PDPN (red), ACTA2 (orange), and HNF4A (blue) are highlighted. j Scatterplot showing the Spearman rho for top CD3E or ACTA2/PDPN- correlated genes to positive relapse recurrence status in the CRC1 and CRC2 subtypes, respectively, in the fresh-frozen FF1–FF2 cohort (left) and FFPE1 (right) cohort. Median correlation is indicated by a black bar. 8 NATURE COMMUNICATIONS | (2020)11:2025 | https://doi.org/10.1038/s41467-020-16000-6 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-16000-6 ARTICLE and EMT pattern was associated with the CIN prognosis sig- number profile, calculated from the methylation array signal nature (Figs. 3e and 4c, g). Furthermore, we compared MCSs for intensity45. This improves cost-effectiveness and makes Meth- TNM stage I (favorable prognosis) to stage IV tumors (poor CORR attractive for analysis of archival FFPE material, where prognosis) in the COREAD and SYSCOL cohorts. Here, the RNA profiling can be difficult6–9. The MethCORR concept bears relative change in MCSs between TNM stages also pointed to a resemblance to transcriptome-wide association studies, where relative loss of immune cells and increase in CAF content in late- gene expression is correlated to genetic variation. However, stage, poor prognosis CRC (Fig. 4h). Collectively, the MethCORR MethCORR allows the expression of many more genes to be analysis of seven published prognostic signatures hereby sug- modeled, which indicates that gene expression is stronger asso- gested that poor prognosis is associated with low T-cell content, ciated with DNA methylation than genetic variation46,47. particularly in the immune-infiltrated CRC1 subtype (Fig. 4f), or The high number of MethCORR genes with inferred expres- high CAF content and inflammation-EMT, particularly in the sion may be surprising, as several previous studies reported more immune-depleted CRC2 subtype (Fig. 4g). To investigate the infrequent correlations, when investigating associations between predictions of prognostic cell types in our FF and FFPE cohorts, gene expression and methylation at local enhancers, promoters, we selected the three biomarkers CD3E, ACTA2, and PDPN. and gene bodies20–22. MethCORR instead performs correlation These are well-known markers for T cells42, CAF/myofibro- analysis genome-wide and hereby identify far more associations blasts43, and inflammation-EMT44, respectively, and their most from which expression information can be inferred. Indeed, closely CpG site-associated genes overlapped with regions high- expression-correlated CpGs were often located far from the gene lighted by the prognostic classifiers (compare Fig. 4b, e, f, g, i; locus, in regions with cell-type-specific methylation (Supple- Supplementary Fig. 6). Indeed, top CD3E-associated genes mentary Fig. 4). Hence, MethCORR benefits from associating negatively correlated with patient recurrence status in the cell-type-specific gene expression with cell-type-specific DNA CRC1 subtype and ACTA2/PDPN-associated genes positively methylation patterns to infer expression information for many correlated to patient recurrence in CRC2 (Fig. 4j). genes, even if associations are not functionally linked. Such indirect associations are expected in heterogenous cancer sam- ples, which vary in their content of cancerous and non-cancerous DNA methylation-based biomarkers for CRC prognostication. cell types2–4,48. Support for a genome-wide correlation strategy is To derive DNA methylation biomarkers for the above prognostic also found in two previous studies, which on a smaller scale, cell types we exploited the cell type-specificity of DNA methy- performed RNA expression-correlation analysis with more dis- lation. Comprehensive comparison of multiple cell types identi- tantly located CpGs49,50. However, these studies only included fied low methylation of CpGs within the CD3E, ACTA2, and ~500 CpG sites distributed across the genome compared with PDPN promoter as biomarkers for T cells, CAFs/myofibroblasts, 480,000 sites utilized in MethCORR, and consequently found and inflammation-EMT, respectively (Fig. 5a; Supplementary much fewer strong correlations. Data 13). Indeed, analysis of promoter CpGs in CRC samples MethCORR introduces an expression-correlated measure, the showed that high methylation of the CD3E promoter, reflecting MCS, which enabled identification of the same two CRC subtypes low levels of T-cell infiltration, associated with significantly in all four cohorts analyzed, and this independent of the analyzed poorer RFS in CRC1 in both FF and FFPE cohorts (Fig. 5b). In tissue being FF or FFPE. The subtypes resemble the two major addition, low ACTA2/PDPN promoter methylation, reflecting carcinogenesis pathways described in CRC32 that are character- high CAF/EMT levels, associated with poor RFS in CRC2 ized by epithelial-cell hyper-methylation or chromosomal (Fig. 5b). The biomarkers were superior predictors of RFS as instability (Figs. 2 and 3). We speculate that MethCORR identi- compared with TNM staging and MSI status (Fig. 5c, Supple- fied these well-established carcinogenesis pathways due to the mentary Fig. 7a, b), and the biomarkers were only prognostic relative emphasis of MCSs on cancer epithelial traits over stroma- within the intended subtype (Supplementary Fig. 7c). Finally, to related traits (Supplementary Fig. 3a, b). Also, we observed higher provide a cost-effective alternative to genome-wide methylome correlations between MCSs profiles for matched FF and FFPE analysis, we evaluated CD3E, ACTA2, and PDPN promoter biopsies taken from the same tumor than between RNA and methylation using quantitative methylation-specific PCR (QMSP) iRNA profiles (Table 1). We therefore speculate that MCS-based assays. In addition, a QMSP assay targeting the HNF4A promoter characterization and subtyping is more independent of sample was included for CRC subtyping; HNF4A is upregulated in CRC2 preservation type, which now require further testing. (Fig. 4i) and correspondingly, its promoter is less methylated MethCORR also introduces a map that visualizes genome- in CRC2 (Fig. 5a). We applied our four biomarker assays to wide associations between gene expression and DNA methy- FFPE1 cohort samples, stratified patients into CRC1 and CRC2 lation in CRC (Fig. 3). We envision that MethCORR map using the HNF4A QMSP assay (Fig. 5d), and used CD3E and analysis may provide a framework for more detailed char- ACTA2/PDPN assays as prognostic biomarkers in CRC1 acterization of FF and archival FFPE samples than categorical and CRC2. RFS analysis confirmed that the QMSP assays allowed subtyping alone, e.g., to reveal cellular sources of inter-tumor subtype-specific prognostication using FFPE samples (Fig. 5e and heterogeneity (Fig. 3). In particular, we illustrated that the Supplementary Fig. 7d). MethCORR map can help identify cell types associated with RNA signatures (Figs. 3 and 4) and hereby help to derive DNA Discussion methylation-based biomarkers suitable for FFPE samples We here introduce MethCORR as an approach for uniform (Fig. 5). Our MethCORR map analysis of several prognostic molecular analysis of FF and FFPE samples based on DNA RNA signatures (Fig. 4) showed that they all predicted cancer methylation profiling. MethCORR allows inference of expression aggressiveness to be associated with cell types within the TME: information from DNA methylation for a large number of genes In particular, a high CAF content, inflammation-associated (>11,000; Fig. 1). The inferred expression profiles support iden- EMT, and low T-cell content were associated with poor prog- tical subtype discovery, characterization, and prognostication in nosis (Fig. 4). This agrees with clinically promising biomarkers FF and FFPE cohorts (Figs. 2–5). Notably, MethCORR allows such as the Immunoscore42 and Tumor-Stroma Ratio51. Our three layers of information to be extracted from a DNA methy- analysis of CRC subtype-specific prognostic RNA signatures lation array experiment, namely an inferred gene expression offered additional resolution: the T-cell content was primarily profile, a DNA methylation profile and a chromosome copy- prognostic within the immune-infiltrated CRC1 subtype, NATURE COMMUNICATIONS | (2020)11:2025 | https://doi.org/10.1038/s41467-020-16000-6 | www.nature.com/naturecommunications 9 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-16000-6 a b FF1-FF2 CRC1 FF1-FF2 CRC2 patients CD3E ACTA2 PDPN HNF4A patients stratified by CD3E stratified by ACTA2/PDPN B-lymphocytes (CD19+) 1.00 1.00 Proportion patients Proportion patients T-lymphocytes (CD3+) w/o recurrence w/o recurrence 0.75 0.75 T-lymphocytes (CD8+) 0.50 0.50 NK cells (CD56+) Monocytes (CD14+) 0.25 P = 0.0372 0.25 P = 0.0198 Neutrophils Methylated HR = 2.49, 95% CI [1.02–6.06] 0.00 0.00 HR = 2.14, 95% CI [1.11–4.11] MSC 1 0 50 100 150 0 50 100 150 Smooth muscle cells -value Time (Months) Time (Months) Endothelial cells 0.5 Number at risk Number at risk CD3E Low meth. 62 37 12 0 ACTA2/PDPN Low meth. 77 40 10 0 CAFs 0 CD3E High meth. 61 29 8 1 ACTA2/PDPN High meth. 87 51 18 0 IMR90 fibroblasts Well diff. Adenomas Unmethylated CD3E Low meth. ACTA2/PDPN Low meth. Medium diff. Adenomas CD3E High meth. ACTA2/PDPN High meth. Poor diff. Adenomas COREAD CRC TNM st. II-III FFPE1 CRC1 patients FFPE1 CRC2 patients MethCORR CRC1 stratified by CD3E stratified by ACTA2/PDPN MethCORR CRC2 1.00 1.00 Proportion patients Proportion patients w/o recurrence w/o recurrence 0.75 0.75 CpG site methylation level 0.50 0.50 0.25 0.25 P = 0.0017 P = 0.0125 HR = 5.43, 95% CI [1.22–24.08] HR = 5.90, 95% CI [1.68–20.76] c FF1–FF2 0.00 0 10 20 30 40 0.00 0 10 20 30 40 Univariate Multivariate Time (Months) Time (Months) P value HR (95% CI) P value HR (95% CI) Number at risk Number at risk MethCORR CD3E Low meth. 21 20 19 19 0 ACTA2/PDPN Low meth. 25 23 16 13 0 0.002 2.32 (1.37–3.93) 0.006 2.12 (1.24–3.61) CD3E High meth. 30 26 19 18 0 ACTA2/PDPN High meth. 28 27 25 25 0 (High vs. Low risk) TNM stage CD3E Low meth. ACTA2/PDPN Low meth. 0.047 1.66 (1.01–2.73) 0.141 1.46 (0.88–2.43) (III vs. II) ACTA2/PDPN High meth. CD3E High meth. MSI 0.078 2.28 (0.91–5.69) 0.112 2.10 (0.84–5.24) (MSS vs. MSI) d 100% e CRC1 patients stratified FFPE1 MethCORR subtype Univariate Multivariate 80% by CD3E (QMSP) percentage P value HR (95% CI) P value HR (95% CI) 60% 1.00 MethCORR Proportion patients 0.000 5.65 (2.17–14.74) 0.001 4.84 (1.84–12.72) w/o recurrence (High vs. Low risk) 40% 0.75 TNM stage 20% 0.50 0.011 2.83 (1.26–6.33) 0.063 2.17 (0.96–4.91) (III vs. II) 0% 0.25 MSI CRC1 predicted CRC2 predicted 0.084 3.54 (0.84–14.85) 0.212 2.52 (0.59–10.76) (MSS vs. MSI) by MethCORR by MethCORR P = 0.0023 0.00 (n = 45) (n = 40) 0 10 20 30 40 Time (Months) FF1-FF2-FFPE1 CRC1 called by CRC2 called by Number at risk Univariate Multivariate QMSP QMSP CD3E Low Ct 25 25 24 24 0 P value HR (95% CI) P value HR (95% CI) CD3E High Ct 16 13 10 10 0 MethCORR 0.000 3.00 (1.91–4.72) 0.000 2.66 (1.67–4.22) CD3E Low Ct (High vs. Low risk) TNM stage CD3E High Ct 0.001 2.01 (1.33–3.05) 0.013 1.71 (1.12–2.60) (III vs. II) MSI CRC2 patients stratified 0.017 2.56 (1.19–5.54) 0.041 2.24 (1.03–4.85) (MSS vs. MSI) by ACTA2/PDPN (QMSP) Proportion patients 1.00 w/o recurrence 0.75 0.50 0.25 0.00 P = 0.0042 0 10 20 30 40 Time (Months) Number at risk ACTA2+PDPN Low Ct 10 10 7 5 0 ACTA2+PDPN High Ct 34 33 31 31 0 ACTA2+PDPN Low Ct ACTA2+PDPN High Ct Fig. 5 Validation of subtype-specific prognostic biomarkers in fresh-frozen and FFPE cohorts. a Dot plot showing the methylation levels (β-values) of a CpG site in the promoter region of CD3E, ACTA2, PDPN, and HNF4A in selected cell types, adenomas and CRC samples as evaluated by the Infinium HumanMethylation450 BeadChip array. High and low methylation levels are indicated in red and blue colors, respectively. See Supplementary Data 13 for details of included cell types such as mesenchymal stromal/stem cells (MSCs), natural killer (NK) cells, and cancer-associated fibroblasts (CAFs). b Kaplan–Meier plot showing the relapse-free survival of patients stratified by the CpG methylation level of the CD3E promoter in CRC1 and by the average CpG methylation level of the ACTA2/PDPN promoter in CRC2 of the combined FF1-FF2- and the FFPE1 cohorts. P values (log-rank test) and HR95% CI are indicated. The same β-value cutoff was used in both cohorts (different cutoff for the subtype-specific biomarkers). c Table showing an uni- and multivariate cox regression analysis with MethCORR high and low relapse risk groups (a high relapse risk group was samples with high CD3E methylation levels in CRC1 or low average ACTA2/PDPN methylation levels in CRC2), TNM stage, and MSI status in the combined FF1–FF2 cohort, the FFPE1 cohort, and all cohorts combined. d Histogram showing the overlap in CRC1 and CRC2 status prediction by NMF clustering using MCSs or by QMSP in FFPE samples from the FFPE1 cohort. e Kaplan–Meier plot showing the relapse-free survival of CRC1 patients stratified by CD3E QMSP assay ΔCt values and in CRC2 by ACTA2/PDPN QMSP ΔCt-values in a total of 85 FFPE samples from the FFPE1 cohort. P values (log-rank test) and HR95% CI are indicated. whereas CAF-content/inflammation-EMT was only prognostic application of the four QMSP assays in CRC samples confirmed in the less immune-infiltrated CRC2 subtype (Fig. 5). This and reproduced the RFS analysis derived from the more costly supported our previous observations of subtype-specific prog- DNA methylome profiles (Fig. 5). Collectively, this illustrates nostic biomarkers3. To aid further testing of subtype-specific the ability of MethCORR to help derive DNA methylation prognostication, we established four simple QMSP assays for biomarkers from transcriptional signatures by extracting cell- cost-efficient CRC subtyping and prognostication. The type information from their expression-correlated CpGs. 10 NATURE COMMUNICATIONS | (2020)11:2025 | https://doi.org/10.1038/s41467-020-16000-6 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-16000-6 ARTICLE Finally, MethCORR can provide high-quality gene expression review boards and ethical committees and written informed consent was obtained measures in samples with poor RNA quality, such as archival from all patients. The FFPE2 cohort (IDIBELL) encompasses 56 samples collected at Medical Oncology Service of ICO Badalona-Germans Trias i Pujol Research FFPE samples for which confident RNA profiling is Institute (IGTP), Spain. None of the patients received neoadjuvant therapy. The challenging6–9. Our analysis of matched FFPE and FF tissue tumors were histologically classified and staged according to the UICC TNM showed that iRNA expression profiles from FFPE tissue resem- staging system. Cancer cell percentage was evaluated individually by two trained bled the RNA-sequencing profiles of the FF tissue better than researchers, and when necessary, tumor biopsies were macroscopically trimmed to RNA-sequencing profiles of the FFPE tissue. In PCA, matched enrich the fraction of neoplastic cells. Patients were followed according to the national clinical guidelines and written informed consent was obtained from all FFPE iRNA and FF RNA-sequencing profiles clustered sample patients. Clinical information regarding the COREAD, SYSCOL, COLOFOL, and wise, while matched RNA-sequencing profiles of FFPE and FF IDIBELL cohort samples is presented in Supplementary Table 1. tissue clustered according to preservation type. Preservation type- dependent clustering of FFPE and FF RNA-sequencing profiles DNA methylome data. FF tumors from the SYSCOL cohort were macrodissected have been reported previously, even in studies that report very to enrich the fraction of neoplastic cells and DNA was extracted from serial high correlation between RNA-sequencing profiles of matched cryosections using the Puregene DNA purification kit (Gentra Systems). Integrity FFPE and FF samples52,53. We acknowledge that recent studies of the genomic DNA from FF samples was assessed by 1.3% agarose gel analysis and only samples containing a high molecular weight smear (~50 KDa) were focusing on newly produced FFPE samples with optimal fixation analyzed further. Bisulfite (BS) conversion of 600 ng DNA of each sample was and short storage time have reported improved correlations performed according to the manufacturer’s recommendations for the Illumina between matched FFPE and FF RNA-sequencing profiles53–55. Infinium Assay (EZ DNA methylation kit. Zymo Research. Cat. No. D5004). Next, However, such samples are not standard in the clinical FFPE DNA methylation profiling was performed using Infinium HumanMethylation450 BeadChip technology (HM-450K; Illumina), as described by the manufacturer. archives. A large study, focusing on clinical FFPE samples, stored FFPE tumors from the COLOFOL FFPE1 cohort were macrodissected to enrich for many years, found that gene expression quantification was the fraction of neoplastic cells, DNA was extracted using the QIAamp DNA FFPE achieved in only 60% of samples and that correlation between Tissue kit (Qiagen) and all samples passed the Infinium FFPE quality control biological replicates was very variable8. (Infinium FFPE QC kit, Illumina). For methylation profiling 500 ng DNA The robustness of MethCORR likely reflects that the Illumina underwent FFPE DNA restoration (Infinium HD FFPE DNA restore kit, Illumina) after BS conversion and profiling was performed using Infinium Infinium HumanMethylation microarray produces highly con- HumanMethylationEPIC BeadChip technology (HM-EPIC; Illumina), as described cordant results in FFPE and FF samples when using DNA by the manufacturer. restoration for FFPE samples (Supplementary Fig. 2h)14–17. FFPE tumors from the IDIBELL FFPE2 cohort were macrodissected to enrich Furthermore, the DNA methylation β-values are calculated as the the fraction of neoplastic cells. DNA was extracted using the QIAamp DNA FFPE Tissue kit (Qiagen) and all samples passed the Infinium FFPE quality control ratio between methylated and unmethylated CpG sites at a given (Infinium FFPE QC kit, Illumina). For methylation profiling 250–500 ng DNA genomic position. Hence, although a genomic region is affected underwent FFPE DNA restoration (Infinium HD FFPE DNA restore kit, Illumina) by degradation, the ratio between the methylated and unmethy- after BS-conversion and profiling was performed using the Infinium lated fragments (i.e., the DNA methylation β-value) would HumanMethylation450 BeadChip technology (HM-450K; Illumina) as described by the manufacturer. For both the SYSCOL, FFPE1, and FFPE2 cohort the expectably be robust. By contrast, RNA profiling is highly affected methylation β-values for each CpG site on the BeadChip were derived using the by RNA degradation26 and the RNA quality obtainable from ChAMP R-package57 using the champ.import and champ.norm functions. FFPE is often compromised6–9. In agreement, tumor samples HM-450K DNA methylation profiles of the COREAD samples were acquired with the lowest correlation between iRNA and measured RNA from the UCSC XENA Public Data Hubs24 [https://xena.ucsc.edu/public-hubs/] and the GDC Data Portal25 [https://portal.gdc.cancer.gov/] as normalized DNA expression had lower RNA quality scores than samples with high methylation β-values. Missing β-values were imputed using the R-package correlations, whereas 450K methylation data quality did not differ Impute58. All DNA methylation measurements were performed once for each (Supplementary Fig. 2g). This suggest that expression profiling of distinct sample. FF samples is influenced by even slight RNA degradation, as reported previously26. RNA-sequencing data. FF tumors from the SYSCOL cohort were macrodissected In conclusion, DNA methylation profiling and MethCORR to enrich the fraction of neoplastic cells and total RNA from serial cryosections analysis enables reliable and robust gene expression estimates to were extracted using the RNeasy Mini Kit (Qiagen). RNA integrity was assessed be obtained from clinical samples with compromised RNA using the Agilent RNA 6000 Nano Kit on an Agilent 2100 Bioanalyzer and >98% of analyzed samples had a RNA integrity number (RIN) > 6. Paired end mRNA quality. Furthermore, MethCORR data can be used to obtain sequencing was performed using 500 ng total RNA for library preparation with the clinically relevant information on tumor subtypes, cellular het- TruSeq RNA Sample Prep Kit v2 and the TruSeq SBS Kit v3 was used for erogeneity, and to develop prognostic biomarkers. Consequently, sequencing aiming for a minimum of 40 Million reads per sample. Sequencing MethCORR represents an effective mean to unlock the unique reads were mapped to the human genome issue HG19 (hg19) using the Tophat2 mapper (Tophat: v2.0.1059) and estimating fragments per kilobase of exon per and extensive resource of FFPE tissues in the pathology archives. million fragments mapped (FPKM) values for Ensembl genes using Cufflink We envision that MethCORR in the future will be established for (Cufflinks: v2.2.1; Gencode v15 annotation w/o Pseudogenes60). many other cancer types. RNA-sequencing profiles for the COREAD samples were acquired from the UCSC XENA Public Data Hubs24 [https://xena.ucsc.edu/public-hubs/] as log2(FPKM + 1) normalized RNA expression values for 20,530 genes and via the Methods GDC Data Portal25 [https://portal.gdc.cancer.gov/] as FPKM normalized RNA CRC patient cohorts. The COREAD cohort encompasses mucosa and UICC TNM expression values for 60,483 transcripts. During comparison of RNA-sequencing stage I–IV CRC samples collected as part of TCGA project. All information data from nine matched FF and FFPE samples, only data originating from the same regarding COREAD samples including processed DNA methylation data, RNA TCGA source center (indicated in Supplementary Data 11) were analyzed. expression data, gene-level copy-number data, and clinical patient information Correlations between RNA sequencing in FF and iRNA expression in FFPE (phenotype) were acquired via the UCSC XENA Public Data Hubs24 [https://xena. samples were analyzed using RNA-sequencing data from TCGA source center 22 ucsc.edu/public-hubs/] and the GDC Data Portal25 [https://portal.gdc.cancer.gov/]. (7 of 9 samples; 2 samples from TCGA source 23), as the GDC MethCORR matrix The SYSCOL and FFPE1 cohorts were acquired from the CRC biobank at the used for iRNA calculation was generated using RNA-sequencing data from samples Department of Molecular Medicine, Aarhus University Hospital, Denmark. primarily originating from TCGA source center 22 (76% of samples). All RNA- SYSCOL samples were collected at hospitals in the central region of Jutland, sequencing measurements were performed once for each distinct sample. Denmark from 1999–20133. The FFPE1 cohort encompasses CRC samples from the prospective study COLOFOL56 collected at hospitals in the central region of Jutland, Denmark. None of the patients received neoadjuvant therapy. The tumors Datasets used for MethCORR development. The MethCORR development were histologically classified and staged according to the UICC TNM staging strategy was independently applied in three CRC datasets of paired RNA expres- system. Cancer cell percentage was evaluated individually by two trained sion and DNA methylation data (Supplementary Data 1, 6, and 8) hereby gen- researchers, and when necessary, tumor biopsies were macroscopically trimmed to erating three different MethCORR matrixes and sets of linear regression models. enrich the fraction of neoplastic cells. The SYSCOL and COLOFOL study was Primarily, MethCORR development was performed using Infinium Human- conducted in accordance with Danish law and is approved by local institutional Methylation450K BeadChip (HM-450K) DNA methylation and RNA-sequencing NATURE COMMUNICATIONS | (2020)11:2025 | https://doi.org/10.1038/s41467-020-16000-6 | www.nature.com/naturecommunications 11 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-16000-6 data from 394 samples of the COAD and READ cohorts (COREAD) of the TCGA Establishment and analysis of a MethCORR map. The MethCORR map for the project, acquired in normalized format via the UCSC XENA Public Data Hubs COREAD cohort was created by clustering MethCORR genes according to their (Supplementary Data 1). The analysis was performed using log2(FPKM + 1) nor- overlap in expression-correlated CpGs using Cytoscape V3.2.062 and the applica- malized RNA expression values for all available 20,530 RNAs and DNA methy- tion EnrichmentMap63 (Jaccard + Overlap filtering cutoff 0.126). Only CpGs with lation β-values for the 396,065 CpGs, where β-values were provided by the XENA negatively expression-correlated CpGs from the MethCORR matrix were used for Public Data Hubs24. This analysis generated the COREAD MethCORR matrix identifying the overlap given that inclusion of all expression-correlated CpGs in a (Supplementary Data 3) that is used for calculation of MCSs throughout the single map would complicate interpretation as genes with opposite expression- manuscript, unless otherwise indicated and modeling metrics is reported in Sup- correlation to DNA methylation would cluster together. Genes with no significant plementary Data 2 and 4. Second, the MethCORR approach was applied to RNA- CpG overlap to other genes are not included in the graphical representation of the sequencing (20,336 RNAs) and HM-450K DNA methylation profiles (485,512 MethCORR map for visual simplicity. For interpretation, the MethCORR map was CpGs) from 314 samples of the SYSCOL cohort3 (Supplementary Data 5–7) with overlain with several data types including external DNA methylation data, tran- the aim to validate the performance of the MethCORR approach in an independent scriptionally defined marker genes, gene sets, and signatures. To visualize these cohort. Third, the MethCORR approach was applied to 405 TCGA COREAD diverse data types using the MethCORR map, four types of scores were established samples using RNA expression (17,611 RNAs, these were selected from the original as follows: dataset of 60,483 transcripts as they overlap with the RNAs included in the UCSC For DNA methylation datasets (450K/EPIC arrays), MCSs were first calculated XENA RNA dataset) and DNA methylation data (395,011 CpGs) acquired via the for all samples and two types of scores were used for map visualization. The NCI GDC25 (Supplementary Data 8). This analysis was performed to investigate difference in median MCS z-scores (Δmedian MCS z-score) was used to visualize the impact of RNA normalization methods on MethCORR performance (modeling differences between subtypes encompassing multiple samples (such as between metrics in Supplementary Data 9 and 10) and to generate a GDC data based MethCORR subtypes, CMS subtypes, CRIS subtypes, MSI vs. MSS tumors etc.) MethCORR matrix that was used for analysis of the TCGA FFPE samples included whereas MCS z-scores were used for visualization of differences between individual in this study, as data from these FFPE samples were also acquired via the GDC samples within a cohort. MCS z-scores were calculated for each gene within each database (Supplementary Data 11). investigated cohort by subtracting the cohort mean from an individual sample MCS and dividing the difference by the cohort standard deviation. E.g. for analysis of inter-tumor heterogeneity, MCS z-scores were calculated for each gene within Identification of RNA expression-correlated CpG sites. The CRC cohort was the whole COREAD FF1 cohort. For analysis of the cellular composition of the divided in two discovery sets (sets 1–2, each encompassing 40% of samples), TME cluster, MCS z-scores were calculated from a collection of cell types with whereas a third set was reserved for independent validation (set 3, 20% of the available 450K analysis downloaded from either Marmal-aid64, Gene Expression samples; Fig. 1a and Supplementary Data 1, 6, and 8). Genome-wide correlations Omnibus (GEO)65, or Array express (see Supplementary Table 4 and (Spearman) between the expression of each of the RNAs (log2(FPKM + 1)) and the Supplementary Data 13 for details of included samples; before calculation of MCS DNA methylation β-value of each CpG site were calculated independently in z-scores across all sample types the median MCSs were calculated for similar discovery sets 1 and 2 using the publicly available R function “cor”. All non- significant correlation pairs were discarded (Spearman’s correlation P value < 0.01). sample types, such as technical replicates). For transcriptionally defined marker genes, gene sets, and signatures, two types The remaining expression-correlated CpGs were ranked by their Spearman’s rho in of scores were used for map visualization depending on the data format. For simple each discovery set and next by their rank sum within discovery sets 1 and 2 to gene sets and RNA signatures, defined by only one gene list (e.g., either up or identify top common expression-correlated CpGs. From these lists of ranked CpGs downregulated RNAs), a correlation to median MCS (cMCS) was calculated for specific for each RNA, we selected up to 100 CpGs whose methylation β-value most each MethCORR gene. The cMCSs were calculated as the average Pearson negatively or positively correlated with its expression resulting in lists of ≤200 RNA correlation between the median MCS of the gene set and the MCS of each expression-correlated CpGs for each RNA (depending on the number of MethCORR gene within the FF1, FF2, FFPE1, and FFPE2 cohorts. For complex expression-correlated CpGs in the ranked lists). To ensure analysis robustness, especially in FFPE samples, we excluded all CpG sites that had a detection P value gene sets/signatures, defined by two gene lists (e.g., of both up and downregulated genes), a correlation to median MCS difference score (ΔcMCS) was instead > 0.05 (ChAMP package57) in ≥5% of samples in either the SYSCOL, FFPE1, or calculated for each MethCORR gene. The ΔcMCSs were calculated by subtracting FFPE2 cohort. Top ranking CpGs for all analyzed genes for the TCGA COREAD the cMCSs for the downregulated gene set from the cMCSs for the upregulated cohort (datasets acquired via the UCSC XENA Public Data Hubs) can be found in gene set (ΔcMCS = cMCSupreg. − cMCSdownreg.) for each gene. For visualization, Supplementary Data 3. MethCORR map gene nodes were colored according to these MCS z-scores, ΔMCS z-scores, cMCS, and ΔcMCS as indicated in the text. For map visualization of Calculation of MethCORR scores. For each sample we used the methylation β- published prognostic signatures, cMCS were calculated for the five general (non values of the top ≤200 RNA expression-correlated CpGs (for each gene) to cal- subtype-specific) signatures (CRC-11336, ColoGuideEx37, Oncotype DX38, culate a MCS for all genes with both positively and negatively expression-correlated ColoPrint39, and Tian et al.40), as they are single lists of RNAs associated with poor CpGs using the formula: prognosis CRC (only recurrence score genes from the Oncotype DX panel were analyzed, whereas treatment genes were excluded). For the CRC subtype-specific ! SSC prognosis and CIN prognosis signatures ΔcMCS were calculated, as they are 1 X ≤100 X ≤100 complex signatures encompassing lists of RNAs with high and low expression in MCS ¼ β value pos: correl: CpG probe þ 1  β value neg: correl: CpG probe : ≤200 1 1 aggressive CRC3. The MCS formula calculates the average methylation value of the expression- NMF-based consensus clustering and SubMap analysis. NMF consensus correlated CpG sites specific for each gene. Unless otherwise indicated, the clustering was performed using the R-package NMF66 with MCSs as input. The COREAD MethCORR matrix encompassing expression-correlated CpGs for number of classes was determined by the first distinctive reduction in the cophe- 11,222 genes (Supplementary Data 3; MethCORR genes) was used for calculation netic score and silhouette consensus score67 and samples were classified according of MCSs throughout the manuscript. The use of the MSC formula above and the to consensus class. The similarity of independent subtype predictions was analyzed MethCORR matrix provided in Supplementary Data 3 allow calculation of MCSs using the Genepattern SubMap module (v328,68) using pairwise comparisons of from DNA methylation β-values of any relevant 450K CRC data set of choice. MCSs and the following settings: num. marker genes = 50, number permutations for Fisher’s statistics = 1000, weighted score type = no, null distribution = each. A false discovery rate (FDR) P value < 0.05 was used as significance cutoff (provided Modeling and inferring of RNA expression from MCSs. We modelled the by the Submap software68). relationship between MCSs and RNA expression for each gene in the discovery samples (set 1 + 2; Fig. 1A) using both simple linear (RNA = B0 + B1 × MCS) and polynomial regression models (RNA = B0 + B1 × MCS + B2 × MCS2…+ Bn × CMS and CRIS subtype classification. CMS classification was performed with the MCSn; n = 2–4). The Caret R-package61 was used to perform modeling by 10 × 10- R-package CMSclassifier using the single sample method and nearest CMS as fold cross validation and we used the average RMSE to select the best model for predicted subtype2. RNA expression or iRNA expression were used as input, as each gene. As performances were highly similar for simple linear and polynomial indicated in the text. CRIS classification was performed using the R-package models for most genes, we only selected polynomial models if a ≥5% relative CRISclassifier provided by Isella et al.4 using RNA expression or iRNA expres- decrease in RMSE values were observed over simple linear models. Model per- sion as input, as indicated in the text. formances were independently validated in validation set 3 (Supplementary Data 2, 7, and 9). Genes with well-performing models (R2 > 0.16 in both the discovery (set 1 + 2) and validation (set 3)) were regarded as MethCORR genes and included in Stroma, CIN, DNA methylation, and ESTIMATE scores. Stroma scores for each the MethCORR matrix (Supplementary Data 3), whereas genes with poorer per- gene (fraction of reads of murine origin) was acquired from Isella et al.48. Genes forming models were excluded. For MethCORR genes we inferred RNA (iRNA) with stroma scores >0.5 were considered stromal genes, whereas genes with stroma expression for each gene in each sample using the MCS as input in the gene- scores <0.1 were considered epithelial cancer genes. For the COREAD cohort, specific linear regression models. Information of the gene-specific models are gene- and sample-specific CIN scores were established from the gene-level copy- provided in Supplementary Data 2, which allow calculation of iRNA profiles from number data (GISTIC2 analysis) available at the UCSC XENA Public Data Hubs24. MCSs for any relevant 450K CRC data set of choice. The gene CIN scores were defined for each gene as the standard deviation of the 12 NATURE COMMUNICATIONS | (2020)11:2025 | https://doi.org/10.1038/s41467-020-16000-6 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-16000-6 ARTICLE GISTIC2 copy-numbers of all samples within the COREAD cohort. The sample good clinical annotation and follow-up (Supplementary Table 1). The inclusion CIN scores were defined as the standard deviation of GISTIC copy-number scores criteria were as follow: A minimum of 2 years follow-up and survival after tumor across all genes within a sample calculated for each sample within the COREAD resection, no local recurrence of the disease, no other cancer within 3 years, and no cohort. For non-COREAD cohort samples (without GISTIC2 data) CIN scores synchronous cancers. RFS was measured from date of surgery to verified first were derived from copy-number data extracted from the HM-450K/EPIC radiologic recurrence (distant) and was censored at the last follow-up or death. The methylome BeadChips using the champ.CNA module of the ChAMP R-package57. following average normalized β-value cutoffs were used for the CD3E, ACTA2, and Here, the sample CIN score was defined as the mean interquartile range of the PDPN CpG probes to stratify patients into high- and low relapse risk groups: β- copy-numbers for all chromosomal segments (seg.mean) covered by at least 25 valueCD3E < 1, average β-valueACTA2, PDPN ≤ 1. The following ΔCT cutoffs were used Illumina probes (num. probes). The sample DNA methylation score for each for the CD3E, ACTA2, and PDPN QMSP biomarker assays to stratify patients into sample was defined as the 40th percentile of DNA methylation β-values of all CpG high- and low relapse risk groups: ΔCTCD3E < 19.5, ΔCTACTA2 < 15.55, ΔCTPDPN < sites common to all four CRC cohorts. ESTIMATE Stroma scores and Immune 13.5. Survival analysis was performed using the Kaplan–Meier method with the scores were calculated using the R-package ESTIMATE69 using default parameters Stata/IC 14.2 (StataCorp) software. Significance was evaluated by log-rank test of and MCSs as input. Household gene status was defined as genes included in the list equality. Cox proportional hazards regression analysis was used to assess the of housekeeping genes70 available at [https://www.tau.ac.il/~elieis/HKG/]. impact of MethCORR risk groups, TNM stage, and MSI status on RFS. The pro- portional hazard assumption was tested by a global test of the Schoenfeld residuals. Gene set enrichment analysis. Pre-ranked GSEA was performed using the GSEA 3.0 tool29 using default settings. Genes were pre-ranked according to their Reporting summary. Further information on research design is available in Spearman correlation of their MCS to CRC1 subtype status and gene sets were the Nature Research Reporting Summary linked to this article. considered significantly up- or downregulated for FDR q values < 0.05 (provided by the GSEA software29). The Molecular Signatures Database (MsigDB) gene set collection v6.1 was used with the addition of custom gene sets (Supplementary Data availability Normalized 450K DNA methylation datasets for the TCGA COREAD cohort used in this Table 3). study are publicly available via the UCSC XENA Public Data Hubs24 [https:// xenabrowser.net/datapages/?dataset=TCGA.COADREAD.sampleMap% Immunohistochemistry. Immunohistochemical stainings of CRC tissue sections 2FHumanMethylation450&host=https%3A%2F%2Ftcga.xenahubs. were acquired from the Human Protein Atlas71 [https://www.proteinatlas.org/]. net&removeHub=http%3A%2F%2F127.0.0.1%3A7222] using the “dataset ID: TCGA. The following antibody and tissue sections were chosen (available from v8.pro- COADREAD.sampleMap/HumanMethylation450” and via the GDC Data Portal25 teinatlas.org): ACTA2 (antibody: CAB013531; Pt. 2001, Pt. 1898, Pt. 2468, Pt. [https://portal.gdc.cancer.gov/repository?facetTab=cases&filters=%7B%22op%22%3A% 3074), PDPN (antibody: HPA007534; Pt. 2001, Pt. 1958, Pt. 1898, Pt. 3264), CD3E 22and%22%2C%22content%22%3A%5B%7B%22op%22%3A%22in%22%2C%22content (antibody: HPA043955; Pt. 4724, Pt. 5005, Pt. 4448, Pt. 5004), HNF4A (antibody: %22%3A%7B%22field%22%3A%22cases.project.program.name%22%2C%22value%22% CAB019417; Pt. 2001, Pt. 2151, Pt. 1958, Pt. 3074). 3A%5B%22TCGA%22%5D%7D%7D%2C%7B%22op%22%3A%22in%22%2C% 22content%22%3A%7B%22field%22%3A%22cases.project.project_id%22%2C%22value Identification of cell-type-specific DNA methylation. Genomic regions with %22%3A%5B%22TCGA-COAD%22%2C%22TCGA-READ%22%5D%7D%7D%2C%7B cell-type-specific DNA methylation was identified by comparing multiple cell types %22op%22%3A%22in%22%2C%22content%22%3A%7B%22field%22%3A%22cases. with available 450K analysis downloaded from either Marmal-aid64, GEO65, or samples.portions.is_ffpe%22%2C%22value%22%3A%5B%22false%22%5D%7D%7D%2C Array express72. Median MCSs were calculated for similar sample types, such as %7B%22op%22%3A%22in%22%2C%22content%22%3A%7B%22field%22%3A%22files. technical replicates prior to analysis (see Supplementary Data 13 for details). For data_type%22%2C%22value%22%3A%5B%22Methylation%20Beta%20Value%22%5D% selection of cell-type-specific methylation markers, only CpG probes with a ≥0.3 7D%7D%2C%7B%22op%22%3A%22in%22%2C%22content%22%3A%7B%22field%22% lower methylation β-value in the intended cell type, as compared with other 3A%22files.platform%22%2C%22value%22%3A%5B%22illumina%20human% relevant cell types, were selected. The following genes/CpG probes were included 20methylation%20450%22%5D%7D%7D%5D%7D] as “datatype=methylation beta here: CD3E/cg24612198, ACTA2/cg09990481, PDPN/cg15563963, HNF4A/ value”, “platform=illumina human methylation 450”, and “case/biospecimen filter: cg06640637. samples portions is FFPE=false” for the TCGA-COAD and TCGA-READ project. Normalized RNA-sequencing data sets for the TCGA COREAD cohort used in this study are publicly available via the UCSC XENA Public Data Hubs24 [https://xenabrowser.net/ Quantitative methylation-specific PCR. QMSP was performed using DNA pri- datapages/?dataset=TCGA.COADREAD.sampleMap%2FHiSeqV2&host=https%3A%2F mers specific for unmethylated CD3E, ACTA2, PDPN, and HNF4A gene promoter %2Ftcga.xenahubs.net&removeHub=http%3A%2F%2F127.0.0.1%3A7222] using the regions (See Supplementary Table 5 for primer sequences). BS conversion was “dataset ID: TCGA.COADREAD.sampleMap/HiSeqV2” and via the GDC Data Portal25 performed with the EZ DNA Methylation-Direct™ Kit (ZYMO research) according [https://portal.gdc.cancer.gov/repository?facetTab=cases&filters=%7B%22op%22%3A% to the manufacturer’s protocol. QMSP was performed using the ViiA™ 7 Real-Time 22and%22%2C%22content%22%3A%5B%7B%22op%22%3A%22in%22%2C%22content PCR system (Applied Biosystems). Biomarker assay reactions were carried out in %22%3A%7B%22field%22%3A%22cases.project.program.name%22%2C%22value%22% triplicate in a final volume of 6 µl and contained 2.5 µl TaqMan® Universal PCR 3A%5B%22TCGA%22%5D%7D%7D%2C%7B%22op%22%3A%22in%22%2C% Master Mix, No AmpErase® UNG (Applied Biosystems), 0.15 µl of 20 pmol/µl 22content%22%3A%7B%22field%22%3A%22cases.project.project_id%22%2C%22value forward and reverse primer, 0.2 µl of 5 pmol/µl hydrolysis probe, 0.125 µl TEMPase %22%3A%5B%22TCGA-COAD%22%2C%22TCGA-READ%22%5D%7D%7D%2C%7B hot start DNA polymerase, 0.4 µl of 12.5 pmol/µl dNTP mix, 0.475 µl H2O, and 2 µl of 2.5 ng/µl BS treated DNA template. AluC4A reference gene reactions were %22op%22%3A%22in%22%2C%22content%22%3A%7B%22field%22%3A%22cases. carried out in triplicate in a final volume of 6 µl and contained 2.5 µl TaqMan® samples.portions.is_ffpe%22%2C%22value%22%3A%5B%22false%22%5D%7D%7D%2C Universal PCR Master Mix, No AmpErase® UNG (Applied Biosystems), 0.2 µl of %7B%22op%22%3A%22in%22%2C%22content%22%3A%7B%22field%22%3A%22files. 25 pmol/µl forward and reverse primer, 0.1 µl of 17 pmol/µl hydrolysis probe, analysis.workflow_type%22%2C%22value%22%3A%5B%22HTSeq%20-%20FPKM%22% 0.125 µl TEMPase hot start DNA polymerase, 0.4 µl of 12.5 pmol/µl dNTP mix, 5D%7D%7D%2C%7B%22op%22%3A%22in%22%2C%22content%22%3A%7B%22field 0.475 µl H2O, and 2 µl of 2.5 ng/µl BS treated DNA template. QMSP reactions were %22%3A%22files.experimental_strategy%22%2C%22value%22%3A%5B%22RNA-Seq% mixed in MicroAmp Optical 384 Well reaction Plates (Applied Biosystem) and run 22%5D%7D%7D%5D%7D] as “Experimental strategy=RNA-Seq”, “Workflow Type= on the ViiA™ 7 Real-Time PCR system (Applied Biosystems) with the following HTSeq = FPKM”, and “case/biospecimen filter=samples portions is FFPE=false” for the PCR program: denaturation at 95 °C for 10 min, followed by 40 cycles at 95 °C for TCGA-COAD and TCGA-READ project. 450K DNA methylation and RNA-sequencing 15 s, and 60 °C for 1 min. The ViiA7TM software (Applied Biosystems) was used for data from TCGA CRC patient with matched FF and FFPE samples are publicly avialable evaluation of the fluorescence signals and the ΔCT was calculated by the use of the via the GDC Data Portal25 [https://portal.gdc.cancer.gov/] using the Database UUID reference gene AluC4A. Subtyping was performed using ΔCTHNF4A as a marker for provided in Supplementary Data 11. The RNA-sequencing data from the SYSCOL CRC2 and the average of ΔCTCD3E and ΔCTACTA2 as a marker for CRC1 (high adenoma/carcinoma samples and the SYSCOL 450K methylome data is deposited at stromal/immune cell infiltration). Samples with a ΔCTHNF4A/(ΔCTCD3E + European Genome-phenome Archive (EGA, [https://www.ebi.ac.uk/ega/]), which is ΔCTACTA2)average ratio <0.85 were defined as CRC2 samples. CT values were hosted by the European Bioinformatics Institute (EBI) and the Centre for Genomic measured three times for each sample (technical triplicates). Regulation (CRG). Study accession numbers are: EGAS00001002376 (RNA sequencing) and EGAS00001004293 (methylomes). The dataset and sample ID’s of the other publicly available DNA methylation datasets used in this study are given in Supplementary Statistical analysis and RFS analysis. Unless otherwise noted, statistical sig- Data 13. All other data supporting the findings of this study are available within the nificance of differences between groups was determined using a non-parametric article, its supplementary information files and from the corresponding author upon WRS test. During Submap analysis28 and pre-ranked GSEA29 a FDR-corrected P value < 0.05 was considered significant and P values were provided by the corre- reasonable request. A reporting summary for this article is available as a Supplementary sponding software. During gene list enrichment analysis an adjusted P value < 0.05 Information file. was considered significant (provided by the Enrichr software30). During eFORGE analysis a Q-value < 0.05 was considered significant (provided by the eFORGE Code availability software). RFS analysis was performed in UICC TNM stage II–III samples with R codes for calculation of MCSs and iRNA profiles are available upon request. NATURE COMMUNICATIONS | (2020)11:2025 | https://doi.org/10.1038/s41467-020-16000-6 | www.nature.com/naturecommunications 13 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-16000-6 Received: 3 March 2019; Accepted: 2 April 2020; 28. Hoshida, Y., Brunet, J. P., Tamayo, P., Golub, T. R. & Mesirov, J. P. Subclass mapping: identifying common subtypes in independent disease data sets. PloS ONE 2, e1195 (2007). 29. Subramanian, A. et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc. Natl Acad. Sci. USA 102, 15545–15550 (2005). References 30. Kuleshov, M. V. et al. Enrichr: a comprehensive gene set enrichment analysis 1. Puppa, G., Sonzogni, A., Colombari, R. & Pelosi, G. TNM staging system of web server 2016 update. Nucleic Acids Res. 44, W90–W97 (2016). colorectal carcinoma: a critical appraisal of challenging issues. Arch. Pathol. 31. Conesa-Zamora, P. et al. Methylome profiling reveals functions and genes Lab. Med. 134, 837–852 (2010). which are differentially methylated in serrated compared to conventional 2. Guinney, J. et al. The consensus molecular subtypes of colorectal cancer. Nat. colorectal carcinoma. Clin. Epigeneti. 7, 101 (2015). Med. 21, 1350–1356 (2015). 32. Leggett, B. & Whitehall, V. Role of the serrated pathway in colorectal cancer 3. Bramsen, J. B. et al. Molecular-subtype-specific biomarkers improve pathogenesis. Gastroenterology 138, 2088–2100 (2010). prediction of prognosis in colorectal cancer. Cell Rep. 19, 1268–1280 (2017). 33. Grun, D. et al. Single-cell messenger RNA sequencing reveals rare intestinal 4. Isella, C. et al. Selective analysis of cancer-cell intrinsic transcriptional traits cell types. Nature 525, 251–255 (2015). defines novel clinically relevant subtypes of colorectal cancer. Nat. Commun. 34. Polyak, K. & Weinberg, R. A. Transitions between epithelial and mesenchymal 8, 15107 (2017). states: acquisition of malignant and stem cell traits. Nat. Rev. Cancer 9, 5. Wang, W. et al. Molecular subtyping of colorectal cancer: Recent progress, 265–273 (2009). new challenges and emerging opportunities. Semin. Cancer Biol. https://doi. 35. Eckstein, M., Rea, M. & Fondufe-Mittendorf, Y. N. Transient and permanent org/10.1016/j.semcancer.2018.05.002 (2018). changes in DNA methylation patterns in inorganic arsenic-mediated 6. Esteve-Codina, A. et al. A comparison of RNA-Seq results from paired epithelial-to-mesenchymal transition. Toxicol. Appl. Pharmacol. 331, 6–17 formalin-fixed paraffin-embedded and fresh-frozen glioblastoma tissue (2017). samples. PloS ONE 12, e0170632 (2017). 36. Nguyen, M. N. et al. CRC-113 gene expression signature for predicting 7. Norton, N. et al. Gene expression, single nucleotide variant and fusion prognosis in patients with colorectal cancer. Oncotarget 6, 31674–31692 transcript discovery in archival material from breast tumors. PloS ONE 8, (2015). e81925 (2013). 37. Agesen, T. H. et al. ColoGuideEx: a robust gene classifier specific for stage II 8. Zhao, Y. et al. Robustness of RNA sequencing on older formalin-fixed colorectal cancer prognosis. Gut 61, 1560–1567 (2012). paraffin-embedded tissue from high-grade ovarian serous adenocarcinomas. 38. Webber, E. M., Lin, J. S. & Evelyn, P. W. Oncotype DX tumor gene expression PloS ONE 14, e0216050 (2019). profiling in stage II colon cancer. Application: prognostic, risk prediction. 9. Jones, W. et al. Deleterious effects of formalin-fixation and delays to fixation PLoS Curr. 2, https://doi.org/10.1371/currents.RRN1177 (2010). on RNA and miRNA-Seq profiles. Sci. Rep. 9, 6980 (2019). 39. Salazar, R. et al. Gene expression signature to improve prognosis prediction of 10. Zhang, P., Lehmann, B. D., Shyr, Y. & Guo, Y. The utilization of formalin stage II and III colorectal cancer. J. Clin. Oncol. 29, 17–24 (2011). fixed-paraffin-embedded specimens in high throughput genomic studies. Int J. 40. Tian, X. et al. Recurrence-associated gene signature optimizes recurrence-free Genom. 2017, 1926304 (2017). survival prediction of colorectal cancer. Mol. Oncol. 11, 1544–1560 (2017). 11. Yakovleva, A. et al. Fit for genomic and proteomic purposes: Sampling the 41. Grugan, K. D. et al. Fibroblast-secreted hepatocyte growth factor plays a fitness of nucleic acid and protein derivatives from formalin fixed paraffin functional role in esophageal squamous cell carcinoma invasion. Proc. Natl embedded tissue. PloS ONE 12, e0181756 (2017). Acad. Sci. USA 107, 11026–11031 (2010). 12. Hosein, A. N. et al. Evaluating the repair of DNA derived from formalin-fixed 42. Kwak, Y. et al. Immunoscore encompassing CD3+ and CD8+ T cell densities paraffin-embedded tissues prior to genomic profiling by SNP-CGH analysis. in distant metastasis is a robust prognostic marker for advanced colorectal Lab. Investig. 93, 701–710 (2013). cancer. Oncotarget 7, 81778–81790 (2016). 13. Chen, L. X., Liu, P. F., Evans, T. C. & Ettwiller, L. M. DNA damage is a 43. Togo, S., Polanska, U. M., Horimoto, Y. & Orimo, A. Carcinoma-associated pervasive cause of sequencing errors, directly confounding variant fibroblasts are a promising therapeutic target. Cancers 5, 149–169 (2013). identification. Science 355, 752–75 (2017). 44. Astarita, J. L., Acton, S. E. & Turley, S. J. Podoplanin: emerging functions in 14. de Ruijter, T. C. et al. Formalin-fixed, paraffin-embedded (FFPE) tissue development, the immune system, and cancer. Front Immunol. 3, 283 (2012). epigenomics using Infinium HumanMethylation450 BeadChip assays. Lab. 45. Feber, A. et al. Using high-density DNA methylation arrays to profile copy Investig. 95, 833–842 (2015). number alterations. Genome Biol. 15, R30 (2014). 15. Ohara, K. et al. Feasibility of methylome analysis using small amounts of 46. Gamazon, E. R. et al. A gene-based association method for mapping traits genomic DNA from formalin-fixed paraffin-embedded tissue. Pathol. Int. 68, using reference transcriptome data. Nat. Genet. 47, 1091–1098 (2015). 633–635 (2018). 47. Zhang, W. et al. Integrative transcriptome imputation reveals tissue-specific 16. Moran, S. et al. Validation of DNA methylation profiling in formalin-fixed and shared biological mechanisms mediating susceptibility to complex traits. paraffin-embedded samples using the Infinium HumanMethylation450 Nat. Commun. 10, 3834 (2019). Microarray. Epigenetics 9, 829–833 (2014). 48. Isella, C. et al. Stromal contribution to the colorectal cancer transcriptome. 17. Moran, S., Arribas, C. & Esteller, M. Validation of a DNA methylation Nat. Genet. 47, 312–319 (2015). microarray for 850,000 CpG sites of the human genome enriched in enhancer 49. Thompson, J. A., Christensen, B. C. & Marsit, C. J. Methylation-to-expression sequences. Epigenomics 8, 389–399 (2016). feature models of breast cancer accurately predict overall survival, distant- 18. Lokk, K. et al. DNA methylome profiling of human tissues identifies global recurrence free survival, and pathologic complete response in multiple and tissue-specific methylation patterns. Genome Biol. 15, r54 (2014). cohorts. Sci. Rep. 8, 5190 (2018). 19. Bormann, F. et al. Cell-of-Origin DNA methylation signatures are maintained 50. Thompson, J. A. & Marsit, C. J. A methylation-to-expression feature model for during colorectal carcinogenesis. Cell Rep. 23, 3407–3418 (2018). generating accurate prognostic risk scores and identifying disease targets in 20. Wagner, J. R. et al. The relationship between DNA methylation, genetic and clear cell kidney cancer. Pac. Symp. Biocomput 22, 509–520 (2017). expression inter-individual variation in untransformed human fibroblasts. 51. van Pelt, G. W. et al. Scoring the tumor-stroma ratio in colon cancer: Genome Biol. 15, R37 (2014). procedure and recommendations. Virchows Arch. 473, 405–412 (2018). 21. Kulis, M. et al. Epigenomic analysis detects widespread gene-body DNA 52. Hedegaard, J. et al. Next-generation sequencing of RNA and DNA isolated hypomethylation in chronic lymphocytic leukemia. Nat. Genet. 44, 1236–1242 from paired fresh-frozen and formalin-fixed paraffin-embedded samples of (2012). human cancer and normal tissue. PloS ONE 9, e98187 (2014). 22. Zhong, H., Kim, S., Zhi, D. & Cui, X. Predicting gene expression using DNA 53. Li, J., Fu, C., Speed, T. P., Wang, W. & Symmans, W. F. Accurate RNA methylation in three human populations. PeerJ 7, e6757 (2019). sequencing from formalin-fixed cancer tissue to represent high-quality 23. Horvath, S. DNA methylation age of human tissues and cell types. Genome transcriptome from frozen tissue. JCO Precis. Oncol. 2018, https://doi.org/ Biol. 14, R115 (2013). 10.1200/PO.17.00091 (2018). 24. Goldman, M., Craft, B., Brooks, A.N., Zhu, J., Haussler, D. The UCSC Xena 54. Graw, S. et al. Robust gene expression and mutation analyses of RNA- Platform for cancer genomics data visualization and interpretation. https:// sequencing of formalin-fixed diagnostic tumor samples. Sci. Rep. 5, 12335 doi.org/10.1101/326470 (2018). (2015). 25. Grossman, R. L. et al. Toward a shared vision for cancer genomic data. N. 55. Li, P., Conley, A., Zhang, H. & Kim, H. L. Whole-transcriptome profiling of Engl. J. Med. 375, 1109–1112, https://doi.org/10.1056/NEJMp1607591 (2016). formalin-fixed, paraffin-embedded renal cell carcinoma by RNA-seq. BMC 26. Gallego Romero, I., Pai, A. A., Tung, J. & Gilad, Y. RNA-seq: impact of RNA Genom. 15, 1087 (2014). degradation on transcript quantification. BMC Biol. 12, 42 (2014). 56. Hansdotter Andersson, P. et al. The COLOFOL trial: study design and 27. Vermeulen, J. et al. Measurable impact of RNA quality on gene expression comparison of the study population with the source cancer population. Clin. results from quantitative PCR. Nucleic Acids Res. 39, e63 (2011). Epidemiol. 8, 15–21 (2016). 14 NATURE COMMUNICATIONS | (2020)11:2025 | https://doi.org/10.1038/s41467-020-16000-6 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-16000-6 ARTICLE 57. Morris, T. J. et al. ChAMP: 450k Chip Analysis Methylation Pipeline. Louis-Hansen’s Foundation (17-2-0457), Dansk Kræftforskningsfond (DKF-2017-26 - Bioinformatics 30, 428–430 (2014). (26)), the Knud and Edith Eriksen’s Memorial Foundation, the Neye Foundation, and the 58. Hastie, T., et al. Imputing missing data for gene expression arrays. Stanford Manufacturer Einar Willumsen’s Memorial Foundation (6000073). The Danish Cancer University Statistics Department Technical report (1999). Biobank is acknowledged for biological material. We thank P. Celis, L. Nielsen, L. 59. Kim, D. et al. TopHat2: accurate alignment of transcriptomes in the presence Kjeldsen, B. Devantie, B. Trolle, S. Moran, D. Garcia, and C. Arribas for their technical of insertions, deletions and gene fusions. Genome Biol. 14, R36 (2013). support. The results published here are in part based upon data generated by the TCGA 60. Trapnell, C. et al. Transcript assembly and quantification by RNA-Seq reveals Research Network [https://cancergenome.nih.gov/]. unannotated transcripts and isoform switching during cell differentiation. Nat. Biotechnol. 28, 511–515 (2010). 61. Kuhn, M. Building predictive models in R using the caret package. J. Stat. Author contributions T.B.M., C.L.A., and J.B.B. designed the experiments. T.B.M., M.H.R., J.S., H.O., S.S.A., Softw. 28, 1–26 (2008). J.G., A.M.C., M.C.M., A.H.M., S.L., E.T.D., M.E., C.L.A., and J.B.B. performed the 62. Shannon, P. et al. Cytoscape: a software environment for integrated models of experiments and included patients. T.B.M., M.H.R., C.L.A., and J.B.B. analyzed and biomolecular interaction networks. Genome Res. 13, 2498–2504 (2003). interpreted the data. T.B.M., C.L.A., and J.B.B. drafted the manuscript. All authors 63. Merico, D., Isserlin, R., Stueker, O., Emili, A. & Bader, G. D. Enrichment map: reviewed and approved the final manuscript. a network-based method for gene-set enrichment visualization and interpretation. PloS ONE 5, e13984 (2010). 64. Lowe, R. & Rakyan, V. K. Marmal-aid–a database for Infinium Competing interests HumanMethylation450. BMC Bioinform. 14, 359 (2013). The authors declare no competing interests. 65. Barrett, T. et al. NCBI GEO: archive for functional genomics data sets–update. Nucleic Acids Res. 41, D991–D995 (2013). 66. Gaujoux, R. & Seoighe, C. A flexible R package for nonnegative matrix Additional information factorization. BMC Bioinform. 11, 367 (2010). Supplementary information is available for this paper at https://doi.org/10.1038/s41467- 67. Brunet, J. P., Tamayo, P., Golub, T. R. & Mesirov, J. P. Metagenes and 020-16000-6. molecular pattern discovery using matrix factorization. Proc. Natl Acad. Sci. USA 101, 4164–4169 (2004). Correspondence and requests for materials should be addressed to C.L.A. or J.B.B. 68. Reich, M. et al. GenePattern 2.0. Nat. Genet. 38, 500–501 (2006). 69. Yoshihara, K. et al. Inferring tumour purity and stromal and immune cell Peer review information Nature Communications thanks the anonymous reviewer(s) for admixture from expression data. Nat. Commun. 4, 2612 (2013). their contribution to the peer review of this work. Peer reviewer reports are available. 70. Eisenberg, E. & Levanon, E. Y. Human housekeeping genes, revisited. Trends Genet. 29, 569–574 (2013). Reprints and permission information is available at http://www.nature.com/reprints 71. Uhlen, M. et al. Proteomics. Tissue-based map of the human proteome. Science 347, 1260419 (2015). Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in 72. Kolesnikov, N. et al. ArrayExpress update–simplifying data submissions. published maps and institutional affiliations. Nucleic Acids Res. 43, D1113–D1116 (2015). 73. Galko, M. J. & Krasnow, M. A. Cellular and genetic analysis of wound healing in Drosophila larvae. Plos Biol. 2, E239 (2004). Open Access This article is licensed under a Creative Commons 74. Anastassiou, D. et al. Human cancer cells express Slug-based epithelial- Attribution 4.0 International License, which permits use, sharing, mesenchymal transition gene expression signature obtained in vivo. BMC adaptation, distribution and reproduction in any medium or format, as long as you give Cancer 11, 529 (2011). appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless Acknowledgements indicated otherwise in a credit line to the material. If material is not included in the This research is supported by grants from the European Commission FP7 project article’s Creative Commons license and your intended use is not permitted by statutory SYSCOL (UE7-SYSCOL-258236), the Novo Nordisk Foundation (NNF16OC0023182), regulation or exceeds the permitted use, you will need to obtain permission directly from the Danish National Advanced Technology Foundation (056-2010-1), the John and the copyright holder. To view a copy of this license, visit http://creativecommons.org/ Birthe Meyer Foundation, the Danish Council for Independent Research (Medical Sci- licenses/by/4.0/. ences) (DFF − 0602-02128B, DFF – 4183-00619, DFF − 7016-00332B), the Danish Council for Strategic Research (1309-00006B), the Danish Cancer Society (R40- A1965_11_S2, R56-A3110-12-S2, R107-A7035, R133-A8520), the National Cancer © The Author(s) 2020 Institute of the National Institutes of Health (R01 CA207467), the Aage and Johanne NATURE COMMUNICATIONS | (2020)11:2025 | https://doi.org/10.1038/s41467-020-16000-6 | www.nature.com/naturecommunications 15