DNA methylation identifies genetically and prognostically distinct subtypes of myelodysplastic syndromes

Myelodysplastic syndromes (MDS) are a group of diverse hematologic malignancies caused by accumulation of somatic driver mutations in clonally expanded hematopoietic stem and progenitor cells (HSPCs).^1,-3  The progeny of these cells demonstrate impaired myeloid differentiation, resulting in peripheral blood cytopenias, progressive bone marrow failure, and potential progression to acute myeloid leukemia (AML). Genetic lesions are identified in the vast majority of patients with MDS; however, no single mutation profile describes the typical MDS patient.^4,5  Instead, a wide variety of mutation profiles appear capable of generating morphologic and clinical features typical of MDS.⁶  Epigenetic marks, such as DNA methylation (DNAm), are critical for mammalian development and are frequently implicated in oncogenesis.^7,8  Somatic mutations in epigenetic regulators are prevalent in more than half of all MDS cases and include genes involved in maintenance of DNAm (eg, TET2 and DNMT3A)^9,10  and histone modification (eg, ASXL1 and EZH2).^11,12  Epigenetic mutations frequently occur early in MDS pathogenesis,^13,14  and aberrant DNAm has been associated with progression to AML.^15,-17  MDS patients without mutations in epigenetic regulators can show alterations in their epigenome, suggesting a convergent pathogenic mechanism.^18,19  Finally, DNA methyltransferase inhibitors are standard-of-care treatments for certain MDS subtypes because their effects can improve hematopoiesis and prolong overall survival (OS), with associated changes in DNAm and gene expression in responding patients.^20,,,,-25  Because DNAm patterns may be influenced by microenvironmental cues and genetic perturbations, we hypothesized that convergent oncogenic states defined by DNAm may represent a useful tool to understand the different pathobiological mechanisms active in MDS.

Prior studies have attributed altered gene expression to DNA promoter hypomethylation in MDS,^25,-27  and more recently, to aberrant methylation at more distal regulatory regions, such as enhancers.^24,28,,-31  Methylation assays that interrogate genomic regions beyond promoters may better inform which regulatory regions possess biologic importance in the development and progression of MDS. Because the functional relevance of methylation in these regions is not well understood, computational methods that agnostically evaluate DNAm datasets may better classify MDS patients without bias. To characterize biologically relevant features of the MDS methylome, we used a targeted bisulfite padlock probe (BSPP) sequencing method that captures specific nonoverlapping regulatory regions, including differentially methylated regions (DMRs) related to genes involved in pluripotency, differentiation, and cancer,^32,,-35  all known promoters for human National Center for Biotechnology Information (NCBI) Reference Sequence genes as well as all microRNA genes, CCCTC–binding factor binding sites, and DNase I hypersensitive regions.^36,37  We employed a computational approach called the OncoGenic Positioning System (Onco-GPS) that uses nonnegative matrix factorization (NMF) to decompose a high dimensional methylation matrix into components that serve as inputs for consensus clustering, resulting in classification of patients by shared methylation signatures.³⁸  There are several advantages to this approach over standard clustering of raw methylation values. Standard clustering may disproportionately weight the influence of cytosine guanine dinucleotide (CpG) dense regions with highly correlated methylation levels such as CpG islands, whereas NMF first summarizes highly correlated parts of the data into components before clustering, giving equal weighting for regions of lower CpG density, such as enhancers, that may have comparable or even higher biologic importance. NMF decomposition also produces quantitative component amplitudes for each patient and each component, such that we can determine whether the underlying components used for clustering are associated with any clinical or genetic variables themselves, as a means of validating the biological information contained in the components used for clustering. With standard clustering of raw DNAm values, there is no mechanism for this type of validation.

Here, we report an epigenomic analysis using the BSPP methylation platform and Onco-GPS computational approach of 141 genetically characterized MDS patient samples.^1,39  We identify distinct methylation states enriched for specific patterns of somatic mutations, cytogenetic abnormalities, and clinical outcomes, including differences in OS. Our findings suggest that DNAm signatures, similar to genetic mutations and gene expression profiles, may refine our ability to classify and predict clinical outcomes in MDS.

A total of 154 patients with treatment-naive MDS were considered for this study. Samples and data from 140 patients were obtained and processed as previously described.^1,39,-41  Samples and data from an additional 14 patients were included from the University of California San Diego Moores Cancer Center. All samples were collected with patient consent under protocols approved by institutional review boards and in accordance with the Declaration of Helsinki. Clinical data to determine the Revised International Prognostic Scoring System (IPSS-R) score were available at the time of sample collection.⁴² 

Genomic DNA samples were isolated from bone marrow aspirate mononuclear cells, and most were genetically characterized in previous studies.^1,39  Samples collected at University of California San Diego were genetically characterized by targeted capture of DNA and sequenced on the Illumina platform. Alignment and variant calling were performed using a custom pipeline, and variant calls were annotated using validated software utilizing databases of known germline and somatic variants (supplemental Methods), followed by confirmation of these variants by manual review.

BSPP libraries were generated and sequenced using a standard 150-bp paired-end read protocol on Illumina HiSeq 2500. Sequencing reads were trimmed, processed, and aligned as previously described.^36,37  CpGs with ≥10× and <500× read coverage were included and summarized into 25-bp tiles by the coverage-weighted average of CpGs within the 25-bp region.⁴³  Filtering and quality control analysis yielded 141 samples covering 246 088 25-bp tiles for our final cohort (supplemental Methods).

NMF was performed on the 3% most variably methylated CpG tiles (n = 7382) to decompose these methylation data into a smaller matrix of NMF components that served as inputs for consensus clustering.^38,44  NMF is a feature extraction algorithm that combines similar attributes (in this case, methylation values at collections of CpG loci) into summary components that represent cohesive properties of the data (in this case, localized methylation differences among patients), thus reducing redundancy while maintaining interpretability. The number of NMF components and consensus clusters was selected based on the cophenetic correlation coefficient peak and visual inspection of the clustering result after evaluating different numbers of clusters (supplemental Methods; visual abstract). Unsupervised and consensus hierarchical clustering were also performed on the 3% most variably methylated CpG tiles for comparison.

Differential methylation analyses were performed on 25-bp CpG tiles with ≥10× and <500× coverage in >60% of samples in each comparison group, with a minimum of 3 samples per group (R package “MethylSig” version 0.4.4). Differentially methylated CpG tiles were defined as those with false discovery rate (FDR) corrected P < .05 and a mean methylation difference >20% between comparison groups unless otherwise specified.

Gene expression datasets (Affymetrix GeneChip Human Genome U133 Plus 2.0 arrays) and metadata were downloaded from NCBI Gene Expression Omnibus with accession number GSE58831. Probe intensity values were normalized using the “gcrma” R package from Bioconductor, and normalized probe intensities were averaged for genes with >1 probe. Comparison groups for survival analysis in the Gerstung et al data were defined based on the average expression of the given genes.⁴⁵ 

Preprocessed level 3 RNA-seq and Illumina 450k methylation array data for the AML cohort (LAML) were downloaded from the The Cancer Genome Atlas (TCGA) data portal. See figure legends of supplemental Figures 13 and 15 for specific analyses performed.

Differences in Kaplan-Meier survival curves were assessed using the log-rank test. Univariate and multivariate models of OS were constructed using Cox proportional hazards regression. Variables with univariate P < .2 were evaluated for inclusion in the final multivariable model, and the final model was constructed using a forward and backward stepwise procedure optimized by the Aikake Information Criterion. Variables with multiple categories (such as IPSS-R categories) were assessed for inclusion in the final multivariable model based on the full variable. An FDR corrected Fisher's exact test was used for all tests of enrichment, including genetic lesions in specific clusters, DMRs within RefSeq genes and chromatin state segments, and gene ontology (GO) terms. GO term enrichment analysis was performed with R package “clusterProfiler” version 3.6.0. Associations between NMF component amplitudes and specific genetic lesions were tested via the information coefficient and an empirical permutation test (n = 100 000 permutations per comparison) to determine statistical significance of the association.⁴⁶ 

Of the 154 patient samples selected for methylation profiling, 10 were removed due to insufficient coverage and 3 were removed after being identified as extreme outliers (supplemental Methods), leaving 141 samples in the final cohort. Demographic and clinical information is provided in Table 1 and supplemental Table 1. The median age was 72 years with 103 (73%) men. Forty-eight patients (34%) were known to subsequently receive hypomethylating agent (HMA) therapy. HMA-treated patients had higher IPSS-R risk, higher bone marrow blast percentages, and lower hemoglobin levels. Demographic characteristics of the final cohort were similar to those in previously published studies from which most of these samples were selected.^1,39 

To identify coherent subgroups defined by DNAm profiles, we employed the Onco-GPS unsupervised computational clustering approach (schema shown in visual abstract). Repeated NMF decompositions with random seeding resulted in the greatest stability with 5 NMF components, indicated by a peak in the cophenetic correlation coefficient (supplemental Figure 1A). Consensus clustering of patients by the amplitudes of these 5 NMF components resulted in stable clustering solutions (as indicated by the cophenetic correlation) for 4 and 5 patient clusters; however, only the 5-cluster solution provided the most distinct clusters as measured by the average silhouette width (supplemental Figure 1C-F). Subdivisions of <5 clusters yielded groups composed of patient samples with heterogeneous component amplitudes (supplemental Figure 1G-I). The 5-cluster solution yielded clusters composed of samples with high amplitudes for a single component, where each cluster was largely defined by a separate component (Figure 1A-B; supplemental Figure 1). For example, patient samples in cluster A were significantly associated with component 1, but not with any other component. For comparison, we also performed consensus clustering of raw methylation data for the most variably methylated CpG tiles and observed that 5 clusters again provided a similarly stable solution (supplemental Figures 2 and 3).

Each of the 5 methylation clusters was significantly enriched for specific genetic and cytogenetic lesions (Figure 1C-D; supplemental 2). Cytogenetic abnormalities were unequally distributed across clusters, with abnormal and complex karyotype patients enriched in clusters A and B and depleted in clusters C through E (cluster A: complex odds ratio [OR] 9.3, FDR <0.1; cluster B: abnormal OR 7.6, FDR <0.01). Alterations of chromosome 7 (OR 9.4, FDR <0.1) and del(5q) (OR 7.0, FDR <0.1) were enriched in cluster B. Total mutation burden was lowest in cluster A, with a mean of 0.8 mutations per patient compared with a mean of 1.5 to 2 mutations per patient for clusters B through E (Figure 1D; supplemental Table 3). Mutations known to convey adverse risk occurred more frequently in clusters B and C. Cluster B had significant enrichment of mutations in TP53 (OR 9.4, FDR <0.05) and U2AF1 (OR 6.4, FDR <0.01). Cluster C showed enrichment of mutations in EZH2 (OR 20.4, FDR <0.01), ASXL1 (OR 5.8, FDR <0.01), and RUNX1 (OR 5.0, FDR <0.05). The splicing factors were mutually exclusive and found in different clusters, with U2AF1 mutations enriched in cluster B (OR 6.4, FDR <0.01), SRSF2 mutations enriched to a less than significant level in cluster C (OR 3.2, FDR >0.1), and SF3B1 mutations enriched in cluster E (OR 4.4, FDR <0.05). TET2 mutated patients were present in several clusters but uniquely enriched in cluster D (OR 3.2, FDR <0.1).

The enrichment of specific mutations within each methylation cluster appeared to be driven by the association of these mutations with the single NMF component that primarily defined each cluster (Figure 2). For example, although samples with mutations in ASXL1 and RUNX1 were enriched in cluster C (defined primarily by NMF component 3), samples assigned to other clusters with mutations in these same genes also possessed higher NMF component 3 scores (Figure 2B). TET2 mutant samples associated strongly with NMF component 4 and were enriched in cluster D, although TET2 mutant samples assigned to clusters A, B, C, and E also had higher component 4 scores than their wild-type counterparts. Using an information coefficient-based permutation test, we found unique associations between NMF components and genetic lesions (Figure 2; supplemental Figures 4 and 5). The most significant component-mutation associations were for component 3-ASXL1 (IC 0.417, P < .001), component 3-RUNX1 (IC 0.351, P < .001), component 3-EZH2 (IC 0.311, P < .001), and component 4-TET2 (IC 0.394, P < .001). The most prevalent splicing factor mutations were each associated with different NMF components (component 3-SRSF2, IC 0.232, P ≤ .01; component 4-SRSF2, IC 0.262, P ≤ .01; component 2-U2AF1, IC 0.257, P ≤ .01; component 5-SF3B1, IC 0.211, P = .02).

Small differences in clinical features were observed between methylation clusters (Figure 3A), including platelet and neutrophil counts. In addition, the median age among cluster A patients was significantly lower at 63 years compared with a median of 72 years for clusters B through E (P = .01). HMA-treated patients were distributed across all clusters with no enrichment by response status (supplemental Figure 6). OS curves for patients in each cluster identified 2 major patterns with those in clusters B and C displaying inferior OS compared with patients in clusters A, D, and E (Figure 3B). We combined clusters with similar median OS into “high” (clusters B and C) and “low” (clusters A, D, and E) cluster risk groups. In univariate analysis, high and low cluster risk groups had statistically significant differences in OS (hazard ratio [HR], 1.95; 95% confidence interval [CI], 1.24-3.08; P < .01) (Table 2; supplemental Table 4). By a multivariate analysis that assessed known prognostic variables in the IPSS-R, DNAm cluster risk group was the most statistically significant single predictor of OS (HR, 2.02; 95% CI 1.25-3.27; P < .01) (Table 2; supplemental Methods). In a second multivariate model that also considered somatic mutations, DNAm cluster risk group remained a significant predictor for OS (HR, 1.6; 95% CI 0.95-2.71 ; P = .08) (Table 2). Notably, methylation cluster risk group was a stronger predictor of OS in IPSS-R lower- vs higher-risk patients and retained prognostic value in HMA-treated patients (supplemental Figures 7 and 8). A subset of patients without known prognostic mutations could also be stratified by methylation cluster risk groups (supplemental Figure 9).^1,2,47 

Although the 5 methylation clusters were defined using a subset of the most variably methylated regions, we identified genome-wide methylation patterns associated with each cluster. Cluster-specific DMRs were extracted by comparing members of each cluster with nonmembers at all CpG tiles with sufficient coverage (Figure 4; supplemental Figure 10). DMRs specific to cluster B and cluster D were almost entirely hypermethylated, whereas DMRs specific to cluster C and cluster E were primarily hypomethylated. Cluster A–specific DMRs were composed of both hypermethylated and hypomethylated CpG tiles. We then examined DMRs associated with specific genetic subgroups of patients (supplemental Figures 10B and 11). Few CpG tiles reached the level of significance for defining DMRs in the TET2 mutant vs nonmutant comparison. This was surprising given the association of TET2 inactivation and hypermethylation reported in previous studies,^{19,24,29,30,48}  and we hypothesized that this may have been due to the genetic heterogeneity in groups simply defined by TET2 mutation presence/absence. To address this heterogeneity, we compared patients who carried only a single mutation with patients without that mutation, and in these comparisons, we observed a hypermethylation phenotype in TET2-only mutated patients, in agreement with the literature. We also wondered whether comparing TET2 mutant vs nonmutants with shared cluster membership could more faithfully recapitulate the differences driven by genetics, because patients with shared cluster membership should share similar global DNAm patterns and benefit from lower within-group variance. Comparison of TET2 mutant vs nonmutants within the same cluster yielded zero DMRs across all clusters that had TET2 mutants (supplemental Figure 10B). However, when comparing TET2-mutated patients from 1 cluster with TET2-mutants from other clusters, much larger methylation differences were seen (supplemental Figure 10B), which were consistent with the global differences seen in cluster-specific DMRs. These results highlight the high heterogeneity in DNAm patterns in TET2-mutant MDS and suggest that patients’ DNAm patterns are more strongly influenced by DNAm cluster membership (and patterns of comutations) as opposed to their status as TET2 mutant vs nonmutant.

To determine if cluster-specific DMRs occurred more frequently at specific types of genomic regulatory regions, we measured overlap between DMRs and regulatory regions defined in the Roadmap Consortium 15-state chromatin state model for human mobilized CD34⁺ cells and calculated enrichment of DMRs in each type of chromatin state (Figure 4; supplemental Figure 12).⁴⁹  Cluster A DMRs were specifically enriched in regions flanking the TSS of genes that are actively transcribed in CD34⁺ cells, with a trend toward hypermethylation. Clusters C and D DMRs were both highly enriched in states predicted to act as enhancers, with opposing trends in methylation. Clusters B and E DMRs were highly enriched for states predicted to be bivalent promoters and bivalent enhancers.

We next performed GO term enrichment analysis of cluster-specific DMGs and identified enrichment for many ontologies related to T-cell activation and differentiation for cluster A DMGs (Figure 5; supplemental Tables 5 and 6). Cluster B was solely enriched for biological processes related to cell-cell adhesion, largely driven by a region of differential methylation centered on the protocadherin γ gene cluster. Cluster D DMGs were enriched for genes involved in neutrophil mediated immunity. The most significant DMG between clusters was Wilms tumor 1 (WT1), which was hypermethylated in cluster B and hypomethylated in cluster E (Figure 5A-B; supplemental Table 6). The majority of DMRs within WT1 were located in a downstream regulatory region that had high levels of H3K27ac, H3K4me1, and H3K4me3 signal in reference HSPC data, which are considered activating histone marks.⁵⁰  Patients in our study who had above average methylation in the 2500-bp region downstream of the WT1 TSS had significantly shorter OS (Figure 5C). We then examined WT1 gene expression and survival data from the Gerstung et al study of MDS patients and found that patients with above average expression of WT1 had significantly shorter OS (Figure 5).⁴⁵  We examined the relationship between methylation and expression of WT1 from an external cohort of 200 patients with AML published by TCGA Research Network⁵¹  and found a positive correlation between hypermethylation and increased gene expression at CpG loci that overlapped with our study, suggesting an association between methylation at these sites and gene expression. The next most significant DMG was CD93, which was hypermethylated at a regulatory region surrounding the TSS in our prognostically lower-risk clusters A, D, and E. We observed a trend for improved survival in patients with above average methylation of this regulatory region, whereas gene expression data from Gerstung et al identified a novel and significant association between survival and CD93 expression (supplemental Figure 14). We observed a significant negative correlation between CD93 expression and DNAm in the TCGA cohort (supplemental Figure 15), implying that adverse risk clusters B and C with lower average methylation may have had higher expression that correlates with the survival results observed in the Gerstung et al data.

In this study, we applied a computational clustering approach based solely on the DNAm profiles of 141 bone marrow DNA samples from patients with MDS to identify subgroups with different DNAm states defined at key regulatory regions of the genome. Using this agnostic approach, 5 methylation subtypes were described that possess unique patterns of enrichment for genetic lesions, differences in regulatory element methylation, and associations with survival that are independent of prognostic clinical variables and somatic mutations. Our results demonstrate that somatic genetics are insufficient to predict the methylation state for a given patient with MDS, and that DNAm profiles may supplement prognostic information in the context of known molecular and clinical variables.

A key conclusion of our study is that in MDS, various mutational profiles can share a common DNAm state even if these states are enriched for specific mutated genes. Somatic mutations alone are not sole determinants of DNAm, which may be integrating additional cues from elements like the microenvironment, differentiation state of cells, or age-related changes. This could reflect preferential selection for mutations in a given preexisting epigenetic state or shared mechanisms engaged by distinct lesions that converge on common epigenetic profiles. For example, we observed that mutually exclusive mutations in EZH2 and SRSF2 were both enriched in cluster C, consistent with studies suggesting these mutations share pathogenic mechanisms.^3,52  We also observed that patients with TET2 mutations were distributed to clusters B, C, and E but enriched only in cluster D. In differential methylation analysis comparing TET2 mutants and nonmutants within a single cluster, we found no DMRs. In contrast, when we compared TET2 mutants within a single cluster to TET2 mutants in other clusters, we observed differences in methylation that were consistent with the global differences between clusters, suggesting that TET2-mutant patients can diverge to different epigenetic states, potentially driven by patterns of coincident mutations or other factors, thus suggesting that DNAm differences cannot be attributed simply to particular mutations or cytogenetic abnormalities. These results highlight the heterogeneity in DNAm patterns among TET2 mutant MDS and show that diverse patterns of somatic mutations can converge toward comparable epigenetic states.

Mutations in U2AF1, SRSF2, and SF3B1 are known to occur in a mutually exclusive manner possibly because they are generally not tolerated as coincident lesions or because they may share a common pathogenic mechanism.^53,54  In our study, mutations of each of these genes were enriched in separate epigenetic clusters and associated with distinct NMF components, indicating that these mutations drive distinct oncogenic states and likely engage unique pathogenic mechanisms.

Another important finding in our study was the enrichment of cluster-specific DMRs located at distal regulatory regions of the hematopoietic stem cell epigenome. Much of the previous work on DNAm in MDS has focused on CpG-rich promoter regions or LINE elements.^48,55,-57  Although we identified DNAm differences at promoters, there were an equal or greater number of differences at nonpromoter regions of the genome. Notably, the TET2-mutated samples enriched in cluster D had the greatest representation of cluster-specific DMRs at enhancer regions as did the loci that defined component 4. These observations are consistent with the putative role of TET2 in maintaining the methylation state at enhancer regions^29,30,58,59  and may help explain why a number of studies in AML and MDS have struggled to find clear relationships between TET2 mutations and DNAm.^{19,24,30,60,-62}  Consideration of the underlying oncogenic state may resolve the relative impact of mutations like those in TET2 whose effects may vary in different contexts. The lack of a relationship between DNAm and TET2 mutations across all clusters in our study may be partially explained by the heterogeneous patterns of comutations present in TET2-mutant patients, which may also strongly impact their global methylation patterns (eg, mutations in ASXL1, RUNX1, or EZH2). Overall, our differential methylation analysis demonstrates that the differences in methylation between genetic subgroups were vastly overshadowed by differences driven by methylation cluster membership.

We did not observe significant differences in cytopenias or bone marrow blast percentages across clusters. However, when examining clinical outcomes, we discovered that DNAm clusters had 2 distinct patterns of prognostic risk. Although these OS differences are modest and require independent validation, cluster risk groups retained their prognostic significance even when accounting for somatic mutations and clinical factors included in the IPSS-R. This suggests that additional prognostic information may be integrated by convergent epigenetic states influenced by, but not exclusive to, somatic mutations. In support of this idea, we showed that even among the subset of patients lacking mutations in genes with known prognostic impact, DNAm cluster risk groups could stratify patients by OS. In addition, we identified several genes with significant differential methylation between clusters, including WT1, which may contribute to the prognostic associations we identified. This is consistent with prior work that has associated WT1 expression with poor prognosis in MDS and AML.^63,,-66  More novel was the association between shorter OS and CD93 hypomethylation in our cohort and overexpression in an external MDS cohort, because CD93 is a transmembrane receptor implicated in the pathogenesis of several malignancies.^67,-69  The role of CD93 in myeloid malignancies is relatively unknown. However, it was shown to be essential to the oncogenic potential of nonquiescent leukemia stem cells in KMT2A-rearranged AML; therefore, this gene may warrant further investigation in MDS.⁷⁰ 

Finally, our findings are consistent with those of Shiozawa et al, which identified 2 major subtypes of MDS with differences in genetic lesions and time to transformation to AML using an analogous unsupervised classification approach on gene-expression data from primary MDS CD34⁺ cells.⁷¹  Although a different genomic analysis was used on a more purified cell population, their study similarly concluded that somatic genetics were insufficient to define clinically relevant disease characteristics in MDS. In the 2 groups defined in their study, mutations in TET2 and SF3B1 were enriched in, but not unique to, the subgroup with longer time to AML transformation, while abnormalities in chromosome 7 and mutations in RUNX1 and TP53, among others, were more frequent in the subgroup with shorter time to AML transformation. These findings are consistent with the genetic abnormalities enriched in the 2 cluster risk groups defined by our study, where mutations in TET2 and SF3B1 were solely enriched in our low-risk clusters D and E, whereas chromosome 7 abnormalities and TP53 and RUNX1 mutations were enriched only in our high-risk clusters B and C.

A potential limitation of our study is that although samples were collected from treatment-naive MDS patients, only a subset of patients were known to have later received HMA therapy. Although HMA-treated patients were distributed across all methylation clusters, there was no signal associated with eventual response. Larger methylation profiling studies of MDS patients with known treatment outcomes are ongoing and may increase our ability to incorporate predictive epigenetic and genetic biomarkers into clinical practice.

In conclusion, we demonstrate that using novel computational methods to agnostically classify primary MDS patients by their DNAm patterns can identify subtypes of MDS characterized by distinct patterns of genetic lesions, regulatory region methylation, and prognostic risk. Our results highlight the importance of DNAm in MDS disease biology as a convergent oncogenic phenotype and support future studies investigating DNAm profiles as clinically relevant biomarkers.

The authors thank Donna Neuberg for her review and suggestions of our prognostic modeling.

This work was supported by National Institutes of Health, National Cancer Institute grants U01-CA217885 (P.T. and H.Y.) and P30-CA023100 (P.T. and H.Y.), National Human Genome Research Institute grant R01-HG009285 (P.T.), and Evans Foundation grants (R.B.). This work was also supported by the American Foundation for Pharmaceutical Education (B.R.).

Contribution: B.R., T.N.T., and R.B. designed the study; T.N.T. and D.D. performed sequencing experiments to generate BSPP libraries; B.R. and D.D. performed sequencing data analysis; B.R. and T.N.T. carried out statistical analyses; D.D. and K.Z. contributed vital reagents for BSPP; H.Y. and P.T. contributed to the design and execution of the Onco-GPS analysis; and B.R., T.N.T., and R.B. wrote the manuscript, which was edited by all authors.

Conflict-of-interest disclosure: R.B. has consulted for Celgene, Astex, AbbVie, Daiichi-Sankyo, Forty Seven, and NeoGenomics; received honoraria from Xian-Janssen and Celgene; and received research funding from Calgene and Takeda. The remaining authors declare no competing financial interests.

Correspondence: Rafael Bejar, University of California San Diego Moores Cancer Center, 3855 Health Sciences Dr, MC 0820, La Jolla, CA 92093-0820; e-mail: rabejar@ucsd.edu.