Expression and methylation patterns partition luminal-A breast tumors into distinct prognostic subgroups
© The Author(s). 2016
Received: 5 November 2015
Accepted: 27 May 2016
Published: 7 July 2016
Breast cancer is a heterogeneous disease comprising several biologically different types, exhibiting diverse responses to treatment. In the past years, gene expression profiling has led to definition of several “intrinsic subtypes” of breast cancer (basal-like, HER2-enriched, luminal-A, luminal-B and normal-like), and microarray based predictors such as PAM50 have been developed. Despite their advantage over traditional histopathological classification, precise identification of breast cancer subtypes, especially within the largest and highly variable luminal-A class, remains a challenge. In this study, we revisited the molecular classification of breast tumors using both expression and methylation data obtained from The Cancer Genome Atlas (TCGA).
Unsupervised clustering was applied on 1148 and 679 breast cancer samples using RNA-Seq and DNA methylation data, respectively. Clusters were evaluated using clinical information and by comparison to PAM50 subtypes. Differentially expressed genes and differentially methylated CpGs were tested for enrichment using various annotation sets. Survival analysis was conducted on the identified clusters using the log-rank test and Cox proportional hazards model.
The clusters in both expression and methylation datasets had only moderate agreement with PAM50 calls, while our partitioning of the luminal samples had better five-year prognostic value than the luminal-A/luminal-B assignment as called by PAM50. Our analysis partitioned the expression profiles of the luminal-A samples into two biologically distinct subgroups exhibiting differential expression of immune-related genes, with one subgroup carrying significantly higher risk for five-year recurrence. Analysis of the luminal-A samples using methylation data identified a cluster of patients with poorer survival, characterized by distinct hyper-methylation of developmental genes. Cox multivariate survival analysis confirmed the prognostic significance of the two partitions after adjustment for commonly used factors such as age and pathological stage.
Modern genomic datasets reveal large heterogeneity among luminal breast tumors. Our analysis of these data provides two prognostic gene sets that dissect and explain tumor variability within the luminal-A subgroup, thus, contributing to the advancement of subtype-specific diagnosis and treatment.
Breast cancer is a heterogeneous disease exhibiting high tumor variability in terms of the underlying biological mechanisms, response to treatment, and overall survival rate . Accurate identification of the unique biological features characterizing each subtype is pivotal for improving our understanding of the disease, identifying subtype-specific biomarkers, targeted drug development, and better prediction of response to treatment.
Originally, therapeutic decisions in breast cancer were guided by clinicopathologic parameters like tumor size, presence of lymph-node/remote metastases, and histological grade. In addition, the status of three immunohistochemistry biomarkers - estrogen receptor (ER), progesterone receptor (PR), and human epidermal growth factor receptor 2 (HER2/ERBB2) allowed the development of targeted therapies and proved predictive of treatment response .
With the emergence of global molecular profiling techniques, large genomic datasets became available for subtype discovery using unsupervised algorithms. By this methodology, breast samples are partitioned into subgroups using clustering algorithms, such as hierarchical clustering  or K-Means, and then subgroup significance is evaluated using the clinical data associated with the samples.
Initially, microarray data were used to define four molecular breast cancer subtypes (basal-like, HER2-enriched, luminal and normal-like) based on characteristic gene expression signatures in correlation with clinical data . These molecular subtypes correlated reasonably well with the immunohistochemical biomarker-based classification. Thus, basal-like samples are mostly triple-negative (ER-/PR-/Her2-), luminal samples are mostly ER+, and HER2-enriched tumors are characterized by amplification and high expression of the HER2/ERBB2 gene [5, 6].
Subsequent analysis conducted on a larger dataset separated the luminal subtype into two distinct subgroups named luminal-A and luminal-B. Luminal-B tumors have higher expression of proliferation genes including Ki-67, and confer worse prognosis [7–9]. Moreover, luminal-B tumors respond better to chemotherapy, while patients with luminal-A cancer gain most benefit from antiestrogen treatment .
As the partitioning of breast tumors into five molecular subtypes has gained acceptance and popularity, several expression-based predictors have been developed. A central predictor is PAM50, which maps a tumor sample to one of the five subtypes based on the gene expression pattern of 50 genes . Though expected to be more robust than traditional classification systems that rely only on a few biomarkers, the separation between luminal-A and luminal-B by the various predictors is not consistent, suggesting that these molecular subtypes may not represent distinct coherent sample groups .
Other attempts to classify breast tumors were based on other profiling technologies such as miRNA arrays [13, 14], copy number variations  or a combination of several different technologies [16, 17]. The various studies have different levels of agreement with the expression-based molecular subtypes, but taken together they strongly indicate the existence of additional, subtler subtypes than the PAM50 subtypes .
Epigenetic modifications such as DNA methylation arrays, which measure the methylation status of thousands of CpG sites across the genome , were also used for breast cancer classification. DNA methylation changes were shown to play a pivotal role in cancer initiation and progression [20, 21]. Particularly, promoter hyper-methylation was associated with silencing of tumor suppressor genes . Several studies associated breast cancer molecular subtypes with specific methylation patterns , while others showed that methylation data may reveal additional complexity not captured at the expression level, possibly identifying finer patient groups of clinical importance .
The large breast cancer dataset developed and provided by The Cancer Genome Atlas project  includes more than a thousand breast tumor samples characterized by various modern high-throughput genomic technologies. This dataset constitutes a significant leap forward compared to the older microarray-based data. mRNA abundance levels are measured in TCGA dataset using the RNA-Seq technology. This technology has increased sensitivity and a higher dynamic range compared to microarrays [20, 21]. DNA-methylation arrays applied on the same samples can help decipher biological tumor variability by epigenetic modifications not manifested at the gene expression level.
The aim of this study was to improve the classification of breast tumors based on the extensive TCGA expression and methylation data that have recently become available. We utilized these datasets to revisit the current classification of breast tumors into biologically distinct subgroups. Our improved and refined classification may contribute to the precision of diagnosis and thus, to more personalized treatment.
Our initial question was whether unsupervised clustering of all TCGA breast samples using the RNA-Seq data would reconstruct the partition defined by PAM50. As the luminal samples had the highest variability in our global clustering, we also asked how the luminal samples would cluster into two groups based on the RNA-Seq data, how the resulting sample groups would compare to the PAM50 partition into luminal-A and luminal-B, and whether that partition would have a clinical advantage over the PAM50 partition of the luminal samples. Looking into the internal structure of the highly variable luminal-A samples, we asked whether this PAM50 group can be further partitioned into finer subgroups with biological distinctness and clinical significance. We then used enrichment analysis to explore the biological mechanisms underlying the new luminal-A subgroups.
We asked similar questions about breast tumor variability at the epigenetic level. We evaluated the methylation-based partition of all breast tumors, all the luminal samples and the highly heterogeneous luminal-A, and compared the resulting partitions to PAM50. To examine the biological characteristics of differentially methylated CpGs (DMCs) separating the new methylation-based luminal-A subgroups, we conducted enrichment analysis. Finally, we performed multivariate COX survival analysis to determine whether the new subgroups have independent prognostic value.
Data acquisition and preprocessing
TCGA data on invasive carcinoma of the breast were downloaded from the UCSC Cancer Browser web site  together with accompanying clinical information. The downloaded RNA-Seq gene expression dataset (Illumina HiSeq platform, gene level RSEM-normalized , log2 transformed) included 1215 samples of which 11 samples from male patients, 8 metastatic samples, and 30 samples of unknown tissue source were filtered out. PAM50 calls (obtained directly from UNC, including PAM50 proliferation scores) were available for 1148 of the filtered samples, and were distributed as follows: 183 basal-like, 78 HER2-enriched, 534 luminal-A, 203 luminal-B and 150 normal-like.
We also downloaded DNA methylation profiles (Illumina Infinium Human Methylation 450K platform, beta values)  containing 872 samples of which 8 male samples, 5 metastatic samples and 19 samples of unknown tissue source were filtered out. We used only 679 tumor samples for which PAM50 calls were available, including 124 basal-like, 42 HER2-enriched, 378 luminal-A and 135 luminal-B samples. Our analysis used only the 107,639 probes of the Infinium-I design type for which a gene symbol was available. This allowed us to bypass the bias of the two probe designs included on the array, to focus on differentially methylated sites that are associated with known genes, and also to reduce the number of analyzed features.
Unsupervised analysis of the tumor samples
Unsupervised analysis of the various sample subsets was executed by clustering the samples based on the 2000 features (genes or CpGs) showing the highest variability over the samples included in each analysis. We used the K-Means clustering algorithm in Matlab (release 2015a) with correlation distance and 100 replicates from which a solution minimizing the sum of point-to-centroid distances was chosen. Due to the high variability among sample subgroups in the breast cancer datasets, reselecting the top variable genes for the analysis of each sample set (and renormalizing accordingly) is crucial to ensure use of the features most relevant to that set. Each feature was independently centered and normalized over the analyzed samples prior to clustering.
Cohort descriptions for the samples used in each analysis are provided in Additional file 1 (Tables S-1A, S-2A, S-3A for the RNA-Seq analyses and Tables S-6A, S-7A and S-8A for the DNA methylation analysis). The TCGA sample Ids included in each analysis are listed in Additional file 2.
Sample cluster enrichment and survival analysis
To evaluate the clinical relevance of the sample clusters obtained in each unsupervised analysis, we used the extensive clinical information available from TCGA for each sample. Enrichment significance of sample clusters for categorical variables (such as the PAM50 subtype or histological type) was calculated using the false discovery rate (FDR)-corrected hypergeometric test. For numeric variables (such as age, percent tumor nuclei, and others) the difference between sample groups was evaluated using the Wilcoxon rank–sum test (Mann–Whitney U test).
Survival and recurrence-free survival curves were plotted using the Kaplan-Meier estimator  and p values for the difference in survival for each group versus all other groups were calculated using the log-rank (Mantel–Haenszel) test [29, 30]. Cox univariate and multivariate survival analyses were conducted using Matlab implementation; p values were corrected using FDR. The analysis and visualization scripts are publicly available as an interactive graphical tool named PROMO .
Analysis of differentially expressed genes and gene enrichment
A list of genes that have the highest differential expression between the two RNA-Seq-based sample groups LumA-R1 and LumA-R2 was generated by applying the Wilcoxon rank–sum test on all dataset genes exhibiting non-zero variance (n = 19,913) after flooring all dataset values to 1 and ceiling to 14. We selected the 1000 genes exhibiting the most significant p values that also have a median difference of at least 0.5 (log2-transformed RSEM expression values). All genes on the list had significantly higher expression in the LumA-R2 sample group (the lowest p value was 8.1e-28).
Gene enrichment tests were performed on these 1000 genes against a background of all genes included in the rank–sum test. The Expander software suite [32, 33] was used to detect significant enrichments for Gene Ontology (GO) , Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways , Wiki-Pathways  and chromosomal location enrichments. GO tests were also performed using the GOrilla tool . The list of 1000 top differentially expressed genes and detailed results of the enrichment analysis are provided in Additional file 3.
Analysis of differentially methylated CpGs, correlation to expression and CpG enrichment
To identify CpGs that are differentially methylated between LumA-M1 and LumA-M3 samples we applied the rank–sum test on all CpGs that survived our preprocessing and also had non-zero variability in the relevant samples (n = 93,880). We then selected the 1000 CpGs that had the highest significance and a minimal median difference of 0.2 (in Beta values). All selected CpGs had significantly higher mean methylation in the LumA-M1 compared to the LumA-M3 group.
To focus on DMCs with genes that had concomitant changes in expression, we calculated Spearman correlation between each CpG and the expression profile of its associated gene based on the Illumina probe-set annotation. The correlation values enabled the identification of 586 DMCs (rank–sum p value <0.01, median difference >0.2) negatively correlated to expression (R < -0.2) and a second smaller group of 212 DMCs positively correlated (R > 0.2) with expression.
We used the array CpG annotations provided by Illumina to calculate enrichment of each one of the three CpG lists (top 1000 DMCs, 586 negatively correlated DMCs and 212 positively correlated DMCs) for features like differentially methylated regions (DMRs), enhancer regions, UCSC RefGene groups and regulatory feature groups. Gene enrichment analysis was performed on the unique genes composing each CpG list, using the Expander and Gorilla tools as described above. Enrichment for InterPro  terms was calculated using the Database for Annotation, Visualization and Integrated Discovery (DAVID) . Enrichment for tumor suppressor genes was calculated by hypergeometric test based on the TSGene  catalog. The lists of differentially methylated CpGs in addition to detailed results of the enrichment analysis are provided in Additional files 4, 5 and 6.
Separation of luminal-A and luminal-B samples is not reconstructed by RNA-Seq unsupervised analysis
The resulting clusters exhibited moderate correspondence with PAM50 labels: most basal-like, normal and HER2-enriched samples fell into three different clusters (numbers 4, 5, and 3, respectively, listed in decreasing levels of homogeneity), whereas the luminal samples exhibited much greater variability. Importantly, most luminal-A sample were split between two different clusters - a homogenous luminal-A cluster (cluster 2), and a cluster composed of a mix of luminal-A and luminal-B samples (cluster 1).
Furthermore, the samples assigned to cluster 2 exhibited a very distinct expression pattern, overexpressing 1184 genes compared to cluster 1 (out of the 1421 differentially expressed genes, see “Methods”). Cluster 1 samples overexpressed only 229 genes compared to cluster 2 (see Additional file 1: Figure S-1E for per-cluster distribution and Additional file 1: Figure S-1F for results of differential gene expression analysis).
According to these results, the variability within the luminal samples is not sufficiently captured by the PAM50 luminal-A and luminal-B subtypes. Specifically, they suggest that luminal-A samples can be further partitioned into finer subgroups, possibly having clinical meaning.
Unsupervised partition of luminal samples predicts survival and recurrence better than PAM50
Survival analysis performed on the two luminal partitions (the PAM50 luminal-A/luminal-B partition, and the two K-Means clusters shown in Fig. 2a) showed that the RNA-Seq-based clustering partition outperforms the luminal-A/luminal-B distinction in terms of both survival and recurrence (5-year survival plots are shown in Fig. 2b; also see Additional file 1: Figure S-2A for overall survival plots). Hence, the signal identified by our unsupervised analysis of the RNA-Seq data translates into a clinically relevant partition of the luminal samples that has better predictive power than the PAM50 luminal-A/luminal-B partition in terms of both survival and recurrence.
Luminal-A samples have two distinct classes exhibiting clinical significance
The most apparent property of the resulting partition was the general overexpression pattern in LumA-R2 samples compared to LumA-R1 samples. Indeed, out of the 2000 genes selected for clustering, 1276 were differentially expressed and 1068 of them were overexpressed in LumA-R2 samples (based on the FDR-corrected rank-sum test). A very similar partition (chi-square, p = 1.1e-40) with a parallel overexpression pattern was identified on a microarray gene expression dataset also available from TCGA for a subset of the luminal-A samples used here (n = 265). This supports the conclusion that the partition and distinct overexpression pattern we observed are not an artifact originating from RNA-Seq measurement technology or from any normalization protocols applied on the dataset (see Additional file 1, section 4).
The main characteristics distinguishing between the luminal-A subgroups, LumA-R1 and LumA-R2
Ductal (p = 2.1e-05)
Lobular (p = 9.7e-12)
Age, years, average
Tumor nuclei percent
Normal cell percent
Comparing the luminal-A partition shown in Fig. 3a to the groups formed when clustering all the luminal samples (Fig. 2a), we note that almost all LumA-R2 samples are contained within cluster 2 (composed of mainly luminal-A samples), whereas most LumA-R1 are contained within cluster 1 (composed of a mix of luminal-A and luminal-B samples) (see the second label bar in Fig. 3a). This suggests that LumA-R1 samples are more similar in their expression profile to luminal-B samples compared with LumA-R2 samples.
Luminal-A subgroups exhibit distinct immune system expression profiles
The most enriched functional categories among the 1000 genes most differentially expressed between LumA-R1 and LumA-R2 samples
Number of genes
Regulation of immune system process
Immune system process
Regulation of leukocyte activation
Regulation of multicellular organismal process
Regulation of response to external
Regulation of biological quality
Positive regulation of cell activation
T cell activation
Regulation of cell proliferation
Cytokine-cytokine receptor interaction
Hematopoietic cell lineage
Cell adhesion molecules (CAMs)
Chemokine signaling pathway
Complement and coagulation cascades
T cell receptor signaling pathway
Natural killer cell mediated cytotoxicity
Pathways in cancer
TCR signaling pathway
B cell receptor signaling pathway
Complement activation, classical pathway
We also observed that the overexpression of immune receptor genes in LumA-R2 samples was accompanied by overexpression of several chemokine genes (CCL5, CCL17, CCL19, and CCL21) and their corresponding receptors (CCR5, CCR4, and CCR7). Topping the list of overexpressed genes in Lum-A-R2 samples (ranked by p value) is the Interleukin-33 (IL-33) gene, which drives T helper 2 (Th2) responses . In summary, LumA-R2 samples exhibit better prognosis based on several clinical parameters while overexpressing a significant number of genes related to the immune system.
Analysis of DNA methylation identifies a luminal subgroup characterized by hyper-methylation and a significantly poorer outcome
The luminal-A tumors proved to be the most heterogeneous in our gene expression analysis. To further identify and characterize clinically meaningful subgroups within the luminal-A group, we explored breast tumor variability on the epigenetic level as well.
Similar results were obtained when we clustered only the 513 luminal A and B samples (Fig. 5b). Here we used the top 2000 variable genes within these samples, to remove the effect of the other two subtypes on the clustering. Importantly, out of the 127 samples comprising the hyper-methylated cluster 1, which was associated with reduced survival (p = 2.6e-05), 76 samples were labeled as luminal-A, a subtype usually associated with good survival. In other words, approximately 20 % of the 378 luminal-A samples (as called by the expression-based PAM50) included in the analysis, could actually be assigned to a higher risk group based on methylation data (see Additional file 1, section 7 for more details).
The three-way partition by methylation levels and its association to differential survival risk also appeared when we repeated the analysis in the group of 378 luminal-A samples, using the top 2000 variable CpGs on these samples (Fig. 5c). The three methylation-based luminal-A clusters were designated LumA-M1, LumA-M2 and LumA-M3. The 84-sample LumA-M1 cluster (comprising approximately 22 % of the luminal-A samples) was associated with significantly reduced 5-year survival (p = 0.0031).
Furthermore, the methylation-based partitioning of the luminal-A samples (LumA-M1/2/3) correlated significantly with the expression-based partitioning (LumA-R1/2, chi-square p = 4.4e-08). The LumA-M2 cluster was enriched for LumA-R1 samples (p = 1.4e-06) and the LumA-M3 cluster was enriched for LumA-R2 samples (p = 1.6e-08), showing that the expression and the methylation-based patterns are related (see lower bar on Fig. 5c). Overall, we identified a poorer outcome subgroup within the luminal-A subtype, which is distinguished by a robust hyper-methylation pattern.
Analysis of differentially methylated CpGs between the LumA-M1 and LumA-M3 subgroups and their correlation to gene expression
To uncover the biological features characterizing the distinct methylation patterns observed in the luminal-A subgroups, we examined the 1000 top DMCs (see “Methods”) between the hyper-methylated LumA-M1 (n = 84) and the hypo-methylated LumA-M3 (n = 171). These two sample subgroups represent the two extremes of the methylation gradient observed in the luminal-A samples. Of note, all 1000 top DMCs (representing 483 genes) were hyper-methylated in the LumA-M1 samples compared to LumA-M3 samples.
Gene enrichment in the three subsets of CpGs exhibiting differential methylation between the LumA-M1 and LumA-M3 subgroups
Negative: R < -0.2
Positive: R > 0.2
1000 CpGs, 483 genes
586 CpGs, 340 genes
212 CpGs, 125 genes
Anatomical structure development
Pattern specification process
Single organism signaling
Multicellular organismal process
Anatomical structure development
Single multicellular organism process
Cellular developmental process
Single organism developmental process
Single organism signaling
Single organism developmental process
Anatomical structure morphogenesis
Anatomical structure development
Cellular developmental process
Single organism developmental process
Regulation of transcription from RNA polymerase II promoter
Anatomical structure morphogenesis
Single multicellular organism process
Tumor suppressor genes (TSGene 2.0)
AHRR, AKR1B1, BMP2, C2orf40, CDH4, CDO1, CDX2, CNTNAP2, CSMD1, DLK1, DSC3, EBF3, EDNRB, FAT4, FOXA2, FOXC1, GALR1, GREM1, GRIN2A, ID4, IRF4, IRX1, LHX4, MAL, MIR124-2, MIR124-3, MIR125B1, MIR129-2, MIR137, MIR9-3, ONECUT1, OPCML, PAX5, PAX6, PCDH8, PHOX2A, PRKCB, PROX1, PTGDR, RASL10B, SFRP1, SFRP2, SHISA3, SLIT2, SOX7, TBX5, UNC5D, ZIC1
AKR1B1, ASCL1, BIN1, BMP4, CCDC67, CDK6, CDO1, EBF3, GSTP1, ID4, IRX1, L3MBTL4, LRRC4, MAP4K1, MME, NTRK3, PCDH10, PDLIM4, PROX1, PTGDR, RUNX3, SCGB3A1, SFRP1, SLC5A8, SLIT2, UBE2QL1, UNC5B, VIM, WT1
AMH, GATA4, HOPX, HOXB13, LHX4, LHX6, MAP4K1, ONECUT1, PAX5, RASAL1, TBX5, TP73, WT1, ZIC1
The databases Gene Ontology, InterPro and Tumor Suppressor Genes 2.0 were used to test the hyper-methylated genes for enrichment. Group 1 is composed of the 1000 top differentially methylated CpGs with a mean difference of at least 0.2. All the CpGs on this list had significant hyper-methylation in the LumA-M1 samples compared to LumA-M3 samples. Group 2 is composed of the 586 CpGs with a differential methylation p value <0.01, a methylation mean difference >0.2 and Spearman-based correlation with expression <0.2. Group 3 is composed of 212 CpGs with a differential methylation p value <0.01, a methylation mean difference >0.2 and Spearman-based correlation with expression >0.2.
As DNA-methylation is known to regulate gene expression and as hyper-methylation of promoters is associated with gene silencing in cancer , we focused on LumA-M1 hyper-methylated CpGs that affect the expression of their corresponding genes. To this end, we used the RNA-Seq-based expression data available from TCGA for the same 378 analyzed samples to generate a second list of CpGs that are both hyper-methylated in LumA-M1 samples (differential methylation p < 0.01, median difference of 0.2) and that have methylation levels inversely correlated to the expression level of their corresponding gene (Spearman correlation R < -0.2). As can be seen in Table 4, the 586 CpGs that passed this filter (corresponding to 340 genes) had significant over-representation of upstream parts of their corresponding genes (UCSC RefGene Group: TSS and 1st exon p < =4.4e-05) and under-representation of gene body (p = 1.43e-16) and 3'UTR (p = 5.83e-04). In terms of the regulatory feature group, these 586 CpGs had over-representation of "Promoter Associated Cell type specific" elements (p = 1.40e-04) accompanied by highly significant under-representation of "Promoter Associated" elements (p = 2.94e-31), suggesting that the observed hyper-methylation pattern involves tissue-specific promoters. Among the 340 under-expressed genes containing the 586 hyper-methylated CpGs, there were several tumor suppressor genes with under-expression that has previously been observed in breast cancer, such as L3MBTL4 , ID4 , RUNX3 [55, 56], PROX1 , SFRP1  and others. Gene-level and CpG-level enrichment for the negative correlations are shown in column 2 of Tables 3 and 4, respectively.
Interestingly, the 212 LumA-M1 hyper-methylated CpGs that were positively correlated with expression (Spearman R > 0.2) had higher enrichment of development-related GO terms compared with negatively correlated CpGs ("pattern specification process" p = 1.07e-13, "embryonic morphogenesis" p = 1.05e-10, "cell fate commitment" p = 5.49e-10). In contrast to the negatively correlated CpGs, they had high over-representation of "gene body" and under-representation of "TSS" regions (UCSC RefGene Group, p = 9.48e-20 and p = 7.28e-14, respectively). For gene and CpG level enrichment for the positive correlations see column 3 in Tables 3 and 4, respectively.
The differential methylation pattern distinguishing LumA-M1 from LumA-M3 samples could therefore be characterized by hundreds of CpGs that are hyper-methylated in the LumA-M1 samples. Distinct subsets of these CpGs correlate negatively and positively with the expression of developmental genes.
Cox survival analysis
Multivariate Cox analysis of luminal-A subgroups for five-year survival and five-year recurrence
LumA-R (1 vs 2)
LumA-M (2, 3 vs 1)
Age (<60 vs > =60 years)
Pathologic stage (I, II vs. III, IV)
Gene expression profiling has become a useful tool for breast cancer classification and for direction of treatment . Although the HER2-enriched and the basal-like subgroups are well-defined and indicative for anti-Her2 and chemotherapy treatment, respectively, the ER-positive luminal subgroup still presents a clinical challenge. In general, all luminal tumors are candidates for anti-hormonal therapy. However, some tumors within this class, often with a more proliferative potential and conferring poorer outcome, are considered for additional therapy. Accordingly, the common classification based on the molecular intrinsic subtypes divides the luminal tumors into the luminal-A tumors, which have a better outcome, and the more proliferative luminal-B tumor subgroups, which have a worse outcome. However, this classification is sub-optimal for clinical decisions because the luminal tumors present a phenotypic and prognostic range rather than an exact partition to either group.
In this study, we applied unsupervised analysis on breast tumor samples using both expression and methylation profiles to reveal new genetic and epigenetic patterns that correlate with a clinical outcome, and compared them to the PAM50 subtypes. Overall, our analyses showed that the separation between luminal-A and luminal-B (as represented by PAM50 labels) is not clear-cut, but rather represents a phenotypic continuum (as previously observed [12, 60, 61]). In fact, each of the gene expression and methylation datasets used in our analysis separately enabled partitioning of the luminal samples into groups with better prognostic value than that of PAM50.
Furthermore, when we focused on the PAM50-designated luminal-A samples only, the RNA-Seq expression profiles split the luminal-A samples into two subgroups (Fig. 3a). The lobular-enriched LumA-R2 sample group, characterized by a distinct gene over-expression pattern, was associated with significantly reduced recurrence risk compared with the more proliferative LumA-R1 subgroup. Interestingly, genes constituting that over-expression pattern were significantly enriched for functions related to the immune system, including the more specific enrichment of chemokines and genes of upstream T cell receptor signaling pathways. We postulate that the significantly elevated mRNA levels of immune related genes in LumA-R2 samples are indicative of increased infiltration levels of immune system cells into these tumors.
Typically, chemokines serve as ligands that by binding to their corresponding receptors, attract immune system cells to the site where they are secreted [62, 63]. LumA-R2 samples over-expressed several chemokines and their corresponding receptors. The simultaneous over-expression of both the chemokine CCL5 (previously found to be highly expressed by breast cancer cells ) and one of its receptors - CCR5 (expressed among others by CD8+ cytotoxic T cells), suggests that tumor cell-derived CCL5 attracts CD8+ cytotoxic T lymphocytes (CTLs) to LumA-R2 tumors. Similarly, the over-expressed chemokines CCL19 and CCL21 may be expressed by the tumor cells, whereas their CCR7 receptor may be expressed by licensed DC or (less typically) by naive and central memory T cells.
In line with this possibility, the over-expressed genes in LumA-R2 samples included genes typical of CTLs (and also natural killer (NK) cells), which may lead to anti-tumor cytotoxic activities exerted by the granzyme (GZMA and GZMB) and perforin pathways (PRF1). Accordingly, over-expression of T cell activation genes was also detected in patients with LumA-R2 tumors. Notably, the over-expressed genes are concentrated at the upstream part of the T cell receptor-signaling pathway (Fig. 4). At this stage, it is not clear why downstream effectors are not enriched in LumA-R2 samples; however, it is of interest to see that the alpha chain of IL-15R was over-expressed in these samples, suggesting that T cell activation processes may indeed come into effect in this subgroup of patients.
How could the over-expression of the immune genes by LumA-R2 samples be related, if at all, to reduced tumor recurrence? It is possible that only LumA-R2 tumors can release chemoattractants that induce the migration of antigen-specific, possibly beneficial, leukocyte subpopulations to the tumor site. Despite recent reports associating tumor infiltrating lymphocytes with a better prognosis [65–67], it is yet to be determined how enhanced immunogenic activity in the LumA-R2 tumors may improve their outcome. Possibly in the future, this LumA-R2 characteristic pattern may direct emerging immune-checkpoint-related therapies .
The role of epigenetic regulation in malignant processes is increasingly recognized. Indeed, our analysis of DNA methylation data partitioned the breast tumor samples into four clusters showing only moderate agreement with the expression-based PAM50 subtypes. In line with previous studies [24, 69], one cluster had a hypo-methylation pattern and corresponded with the PAM50 basal-like subgroup that was associated with poorer outcome. However, the luminal samples did not cluster neatly into the PAM50 luminal-A and luminal-B subgroups. Instead, three luminal clusters with increasing methylation levels were obtained (clusters 1–3 in Fig. 5a), of which the most hyper-methylated cluster was associated with significantly poorer 5-year prognosis. In fact, even when we clustered only the luminal-A samples (Fig. 5c), the hyper-methylated cluster 1 (LumA-M1) was still associated with significantly poorer survival compared to the other two clusters (LumA-M2 and LumA-M3).
Notably, the top 1000 differentially methylated CpG loci, all hyper-methylated on LumA-M1 samples, had enrichment for genes involved in morphogenesis, differentiation, and developmental processes. Moreover, the CpG hyper-methylation correlated with under-expression of developmental genes, including various tumor suppressor genes. Indeed, hyper-methylation of developmental genes in luminal breast tumors was previously reported [70, 71], secondary to repressive histone marks, which direct de novo methylation. Moreover, hyper-methylation was implicated in normal processes of cell aging and in tumorigenesis . Taken together, the methylation-based analysis suggests a poorer outcome for luminal tumors with a characteristic hyper-methylation pattern, whether in the luminal-A or in the luminal-B subgroups. The hyper-methylation-associated silencing of developmental and tumor suppressor genes may indeed explain these findings. More importantly, within the luminal-A subgroup that is generally associated with a better outcome, the hyper-methylation pattern of the LumA-M1 subgroup marks 84 samples (comprising 22 % of the 378 luminal-A samples) as a high-risk patient group that might benefit from more aggressive treatment.
Last, we showed that the sample partitions induced by the gene expression and DNA methylation patterns are related (p = 4.4e-08; see lower bar in Fig. 5c), mainly because the LumA-M3 samples that are associated with a better outcome are enriched for LumA-R2. However, our attempts to partition the luminal-A samples based on both patterns together did not yield a partition that is better than the separate partitions, in terms of survival prediction or clustering stability. This observation was confirmed by Cox multivariate analysis showing the independent prognostic contribution of each pattern to outcome prediction (Table 5), suggesting that gene expression and methylation hold complementary information, reflecting different aspects of the biological complexity of breast tumors.
Very recently, several novel partitions of luminal breast tumors were proposed [19, 65, 72]. The partitions identified in this study are reinforced by partial though significant similarity to some newly defined groups. LumA-R1 and LumA-R2 clusters are enriched for the proliferative (p = 8.1e-04) and reactive-like (2.4-e04) classes of invasive lobular carcinoma (ILC), respectively, as defined in  (see Additional file 1, section 12). Furthermore, the LumA-M1 cluster is enriched (p = 1.6e-07) for the Epi-LumB group of tumors that are associated with poorer outcome, described by Stefansson et al.  (named Epi-LumB, as it was largely composed of Luminal-B samples, see Additional file 1, section 13). Additional research is needed in order to consolidate the different partitions identified using different procedures into robust and meaningful categories for prognostic and diagnostic use in clinics.
This study emphasizes the large heterogeneity of luminal breast tumors in general, and of luminal-A samples in particular, the inner variability of which was found to be inadequately captured by PAM50 molecular subtypes. Analysis of the RNA-Seq data revealed a partition of the luminal-A samples into groups associated with different risks of 5-year recurrence. We suggest that the over-expression of immune genes in the LumA-R2 group can be ascribed to a higher tendency of its samples to attract tumor-infiltrating lymphocytes, but this requires further research into the mechanism by which the higher infiltrates affect recurrence risk. In the DNA methylation data, a hyper-methylation pattern enriched for developmental genes defined a luminal-A subgroup that was associated with poorer patient survival. In practice, the two prognostic patterns and the lists of genomic features characterizing each of them, can uncover the biological aspects underlying the heterogeneity of luminal-A tumors, improve our ability to classify these tumors into more accurate clinical subgroups, and contribute to the development of novel directed therapies.
CTL, cytotoxic T lymphocytes; DMC, differentially methylated CpG; DMR, differentially methylated region; ER, estrogen receptor; FDR, false discovery rate; GO, Gene Ontology; HER2, human epidermal growth factor 2; IL, interleukin; KEGG, Kyoto Encyclopedia of Genes and Genomes; NK, natural killer; PR, progesterone receptor; TCGA, The Cancer Genome Atlas; TCR, T cell receptor; TNF, tumor necrosis factor
This study was supported in part by the Israeli Science Foundation (grant 317/13), by an IDEA grant of the Dotan Center in Hemato-Oncology, and by the Israeli Center of Research Excellence (I-CORE), Gene Regulation in Complex Human Disease, Center No. 41/11. DN was supported in part by a fellowship from the Edmond J. Safra Center for Bioinformatics, Tel Aviv University.
The authors declare that they have no competing interests.
DN and RS designed the study and analyzed the results. DN conduced the computational analyses. AA, EE, and ABB contributed to the study design and interpreted the results. All authors contributed to writing the manuscript. All authors have read and approved the final version of the manuscript.
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give 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 Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
- Bertos NR, Park M. Review series breast cancer – one term, many entities? J Clin Invest. 2011;121:3789–96.View ArticlePubMedPubMed CentralGoogle Scholar
- Goldhirsch A, Ingle JN, Gelber RD, Coates AS, Thürlimann B, Senn H-J. Thresholds for therapies: highlights of the St Gallen International Expert Consensus on the primary therapy of early breast cancer 2009. Ann Oncol. 2009;20:1319–29.View ArticlePubMedPubMed CentralGoogle Scholar
- Eisen MB, Spellman PT, Brown PO, Botstein D. Cluster analysis and display of genome-wide expression patterns. Proc Natl Acad Sci. 1998;95:14863–8.View ArticlePubMedPubMed CentralGoogle Scholar
- Perou CM, Sørlie T, Eisen MB, van de Rijn M, Jeffrey SS, Rees CA, et al. Molecular portraits of human breast tumours. Nature. 2000;406:747–52.View ArticlePubMedGoogle Scholar
- Reis-Filho JS, Pusztai L. Gene expression profiling in breast cancer: classification, prognostication, and prediction. Lancet. 2011;378:1812–23.View ArticlePubMedGoogle Scholar
- Sotiriou C, Pusztai L. Gene-expression signatures in breast cancer. N Engl J Med. 2009;360:790–800.View ArticlePubMedGoogle Scholar
- Sørlie T, Perou CM, Tibshirani R, Aas T, Geisler S, Johnsen H, et al. Gene expression patterns of breast carcinomas distinguish tumor subclasses with clinical implications. Proc Natl Acad Sci USA. 2001;98:10869–74.View ArticlePubMedPubMed CentralGoogle Scholar
- Sorlie T, Tibshirani R, Parker J, Hastie T, Marron JS, Nobel A, et al. Repeated observation of breast tumor subtypes in independent gene expression data sets. Proc Natl Acad Sci. 2003;100:8418–23.View ArticlePubMedPubMed CentralGoogle Scholar
- Sørlie T. Molecular portraits of breast cancer: tumour subtypes as distinct disease entities. Eur J Cancer. 2004;40:2667–75.View ArticlePubMedGoogle Scholar
- Goldhirsch A, Wood WC, Coates AS, Gelber RD, Thürlimann B, Senn H-J, et al. Strategies for subtypes–dealing with the diversity of breast cancer: highlights of the St. Gallen International Expert Consensus on the Primary Therapy of Early Breast Cancer 2011. Ann Oncol. 2011;22:1736–47.View ArticlePubMedPubMed CentralGoogle Scholar
- Parker JS, Mullins M, Cheang MCU, Leung S, Voduc D, Vickery T, et al. Supervised risk predictor of breast cancer based on intrinsic subtypes. J Clin Oncol. 2009;27:1160–7.View ArticlePubMedPubMed CentralGoogle Scholar
- Weigelt B, Mackay A, A’hern R, Natrajan R, Tan DS, Dowsett M, et al. Breast cancer molecular profiling with single sample predictors: a retrospective analysis. Lancet Oncol. 2010;11:339–49.View ArticlePubMedGoogle Scholar
- Bhattacharyya M, Nath J, Bandyopadhyay S. MicroRNA signatures highlight new breast cancer subtypes. Gene. 2015;556:192–8.View ArticlePubMedGoogle Scholar
- Blenkiron C, Goldstein LD, Thorne NP, Spiteri I, Chin S-F, Dunning MJ, et al. MicroRNA expression profiling of human breast cancer identifies new markers of tumor subtype. Genome Biol. 2007;8:R214.View ArticlePubMedPubMed CentralGoogle Scholar
- Andre F, Job B, Dessen P, Tordai A, Michiels S, Liedtke C, et al. Molecular characterization of breast cancer with high-resolution oligonucleotide comparative genomic hybridization array. Clin Cancer Res. 2009;15:441–51.View ArticlePubMedGoogle Scholar
- Shen R, Olshen AB, Ladanyi M. Integrative clustering of multiple genomic data types using a joint latent variable model with application to breast and lung cancer subtype analysis. Bioinformatics. 2009;25:2906–12.View ArticlePubMedPubMed CentralGoogle Scholar
- Curtis C, Shah SP, Chin S-F, Turashvili G, Rueda OM, Dunning MJ, et al. The genomic and transcriptomic architecture of 2,000 breast tumours reveals novel subgroups. Nature. 2012;486:346–52.PubMedPubMed CentralGoogle Scholar
- Dawson S-J, Rueda OM, Aparicio S, Caldas C. A new genome-driven integrated classification of breast cancer and its implications. EMBO J. 2013;32:617–28.View ArticlePubMedPubMed CentralGoogle Scholar
- Bibikova M, Barnes B, Tsan C, Ho V, Klotzle B, Le JM, et al. High density DNA methylation array with single CpG site resolution. Genomics. 2011;98:288–95.View ArticlePubMedGoogle Scholar
- Jones PA, Baylin SB. The fundamental role of epigenetic events in cancer. Nat Rev Genet. 2002;3:415–28.View ArticlePubMedGoogle Scholar
- Esteller M. Epigenetics in Cancer. N Engl J Med. 2008;358:1148–59.View ArticlePubMedGoogle Scholar
- Baylin SB. DNA methylation and gene silencing in cancer. Nat Clin Pract Oncol. 2005;2 Suppl 1:S4–11.View ArticlePubMedGoogle Scholar
- Holm K, Hegardt C, Staaf J, Vallon-Christersson J, Jönsson G, Olsson H, et al. Molecular subtypes of breast cancer are associated with characteristic DNA methylation patterns. Breast Cancer Res. 2010;12:R36.View ArticlePubMedPubMed CentralGoogle Scholar
- Rønneberg JA, Fleischer T, Solvang HK, Nordgard SH, Edvardsen H, Potapenko I, et al. Methylation profiling with a panel of cancer related genes: association with estrogen receptor, TP53 mutation status and expression subtypes in sporadic breast cancer. Mol Oncol. 2011;5:61–76.View ArticlePubMedGoogle Scholar
- Koboldt DC, Fulton RS, McLellan MD, Schmidt H, Kalicki-Veizer J, McMichael JF, et al. Comprehensive molecular portraits of human breast tumours. Nature. 2012;490:61–70.View ArticleGoogle Scholar
- Goldman M, Craft B, Swatloski T, Cline M, Morozova O, Diekhans M, et al. The UCSC Cancer Genomics Browser: update 2015. Nucleic Acids Res. 2015;43:D812–7. [Internet] [cited 2015 Mar 31]. Available from: http://nar.oxfordjournals.org/content/early/2014/11/11/nar.gku1073.short.
- Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011;12:323.View ArticlePubMedPubMed CentralGoogle Scholar
- Kaplan EL, Meier P. Nonparametric estimation from incomplete observations. J Am Stat Assoc. 1958;53:457–81.View ArticleGoogle Scholar
- Horwitz RI. Statistical aspects of the analysis of data from retrospective studies of disease. J Chronic Dis. 1979;32:ii.Google Scholar
- Bland JM, Altman DG. The logrank test. BMJ. 2004;328:1073.View ArticlePubMedPubMed CentralGoogle Scholar
- Netanely D, Laufer I, Shamir R. PROMO: profiler of multi-omics data [Internet]. Available from: http://acgt.cs.tau.ac.il/promo/
- Shamir R, Maron-Katz A, Tanay A, Linhart C, Steinfeld I, Sharan R, et al. EXPANDER–an integrative program suite for microarray data analysis. BMC Bioinformatics. 2005;6:232.View ArticlePubMedPubMed CentralGoogle Scholar
- Ulitsky I, Maron-Katz A, Shavit S, Sagir D, Linhart C, Elkon R, et al. Expander: from expression microarrays to networks and functions. Nat Protoc. 2010;5:303–22.View ArticlePubMedGoogle Scholar
- Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, et al. Gene Ontology: tool for the unification of biology. Nat Genet. 2000;25:25–9.View ArticlePubMedPubMed CentralGoogle Scholar
- Kanehisa M, Goto S. KEGG: Kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28:27–30.View ArticlePubMedPubMed CentralGoogle Scholar
- Kelder T, van Iersel MP, Hanspers K, Kutmon M, Conklin BR, Evelo CT, et al. WikiPathways: building research communities on biological pathways. Nucleic Acids Res. 2012;40:D1301–7.View ArticlePubMedGoogle Scholar
- Eden E, Navon R, Steinfeld I, Lipson D, Yakhini Z. GOrilla: a tool for discovery and visualization of enriched GO terms in ranked gene lists. BMC Bioinformatics. 2009;10:48.View ArticlePubMedPubMed CentralGoogle Scholar
- Mitchell A, Chang H-Y, Daugherty L, Fraser M, Hunter S, Lopez R, et al. The InterPro protein families database: the classification resource after 15 years. Nucleic Acids Res. 2015;43:D213–21.View ArticlePubMedGoogle Scholar
- Huang DW, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009;4:44–57.View ArticleGoogle Scholar
- Zhao M, Sun J, Zhao Z. TSGene: a web resource for tumor suppressor genes.Nucleic Acids Research. 2013;41(Database issue):D970-D976. doi: 10.1093/nar/gks937.
- Zola H, Swart B, Banham A, Barry S, Beare A, Bensussan A, et al. CD molecules 2006–human cell differentiation molecules. J Immunol Methods. 2007;319:1–5.View ArticlePubMedGoogle Scholar
- Penninger JM, Irie-Sasaki J, Sasaki T, Oliveira-dos-Santos AJ. CD45: new jobs for an old acquaintance. Nat Immunol. 2001;2:389–96.PubMedGoogle Scholar
- Kuhns MS, Davis MM, Garcia KC. Deconstructing the form and function of the TCR/CD3 complex. Immunity. 2006;24:133–9.View ArticlePubMedGoogle Scholar
- Kehrl JH, Riva A, Wilson GL, Thévenin C. Molecular mechanisms regulating CD19, CD20 and CD22 gene expression. Immunol Today. 1994;15:432–6.View ArticlePubMedGoogle Scholar
- Chen L, Flies DB. Molecular mechanisms of T cell co-stimulation and co-inhibition. Nat Rev Immunol. 2013;13:227–42.View ArticlePubMedPubMed CentralGoogle Scholar
- Voskoboinik I, Whisstock JC, Trapani JA. Perforin and granzymes: function, dysfunction and human pathology. Nat Rev Immunol. 2015;15:388–400.View ArticlePubMedGoogle Scholar
- Cronin SJF, Penninger JM. From T-cell activation signals to signaling control of anti-cancer immunity. Immunol Rev. 2007;220:151–68.View ArticlePubMedGoogle Scholar
- Jovanovic IP, Pejnovic NN, Radosavljevic GD, Pantic JM, Milovanovic MZ, Arsenijevic NN, et al. Interleukin-33/ST2 axis promotes breast cancer growth and metastases by facilitating intratumoral accumulation of immunosuppressive and innate lymphoid cells. Int J Cancer. 2014;134:1669–82.View ArticlePubMedGoogle Scholar
- Tommasi S, Karm DL, Wu X, Yen Y, Pfeifer GP. Methylation of homeobox genes is a frequent and early epigenetic event in breast cancer. Breast Cancer Res. 2009;11:R14.View ArticlePubMedPubMed CentralGoogle Scholar
- Abate-Shen C. Deregulated homeobox gene expression in cancer: cause or consequence? Nat Rev Cancer. 2002;2:777–85.View ArticlePubMedGoogle Scholar
- Shah N, Sukumar S. The Hox genes and their roles in oncogenesis. Nat Rev Cancer. 2010;10:361–71.View ArticlePubMedGoogle Scholar
- Herman JG, Baylin SB. Gene silencing in cancer in association with promoter hypermethylation. N Engl J Med. 2003;349:2042–54.View ArticlePubMedGoogle Scholar
- Addou-Klouche L, Adélaïde J, Finetti P, Cervera N, Ferrari A, Bekhouche I, et al. Loss, mutation and deregulation of L3MBTL4 in breast cancers. Mol Cancer. 2010;9:213.View ArticlePubMedPubMed CentralGoogle Scholar
- Noetzel E, Veeck J, Niederacher D, Galm O, Horn F, Hartmann A, et al. Promoter methylation-associated loss of ID4 expression is a marker of tumour recurrence in human breast cancer. BMC Cancer. 2008;8:154.View ArticlePubMedPubMed CentralGoogle Scholar
- Chen L-F. Tumor suppressor function of RUNX3 in breast cancer. J Cell Biochem. 2012;113:1470–7.PubMedPubMed CentralGoogle Scholar
- Huang B, Qu Z, Ong CW, Tsang Y-HN, Xiao G, Shapiro D, et al. RUNX3 acts as a tumor suppressor in breast cancer by targeting estrogen receptor α. Oncogene. 2012;31:527–34.View ArticlePubMedGoogle Scholar
- Versmold B, Felsberg J, Mikeska T, Ehrentraut D, Kohler J, Hampl JA, et al. Epigenetic silencing of the candidate tumor suppressor gene PROX1 in sporadic breast cancer. Int J Cancer. 2007;121:547–54.View ArticlePubMedGoogle Scholar
- Klopocki E, Kristiansen G, Wild P, Klaman I, Castanos-Velez E, Singer G, et al. Loss of SFRP1 is associated with breast cancer progression and poor prognosis in early stage tumors. Int J Oncol. 2004;25:641–9.PubMedGoogle Scholar
- Prat A, Pineda E, Adamo B, Galván P, Fernández A, et al. Clinical implications of the intrinsic molecular subtypes of breast cancer. Breast. 2015;24:S26–35.View ArticlePubMedGoogle Scholar
- Alizart M, Saunus J, Cummings M, Lakhani SR. Molecular classification of breast carcinoma. Diagnostic Histopathol. 2012;18:97–103.View ArticleGoogle Scholar
- Weigelt B, Pusztai L, Ashworth A, Reis-Filho JS. Challenges translating breast cancer gene signatures into the clinic. Nat Rev Clin Oncol. 2012;9:58–64.View ArticleGoogle Scholar
- Balkwill F. Cancer and the chemokine network. Nat Rev Cancer. 2004;4:540–50.View ArticlePubMedGoogle Scholar
- Balkwill F. Chemokine biology in cancer. Semin Immunol. 2003;15:49–55.View ArticlePubMedGoogle Scholar
- Luboshits G, Shina S, Kaplan O, Chaitchik S, Keydar I, Ben-baruch A. Elevated expression of the CC chemokine regulated on activation, normal t cell expressed and secreted (RANTES) in advanced breast carcinoma. Cancer Res. 1999;59:4681–7.PubMedGoogle Scholar
- Salgado R, Denkert C, Demaria S, Sirtaine N, Klauschen F, Pruneri G, et al. The evaluation of tumor-infiltrating lymphocytes (TILs) in breast cancer: recommendations by an International TILs Working Group 2014. Ann Oncol. 2014;26:259–71.View ArticlePubMedGoogle Scholar
- Denkert C, Loibl S, Noske A, Roller M, Müller BM, Komor M, et al. Tumor-associated lymphocytes as an independent predictor of response to neoadjuvant chemotherapy in breast cancer. J Clin Oncol. 2010;28:105–13.View ArticlePubMedGoogle Scholar
- Denkert C. The immunogenicity of breast cancer--molecular subtypes matter. Ann Oncol. 2014;25:1453–5.View ArticlePubMedGoogle Scholar
- Bedognetti D, Hendrickx W, Marincola FM, Miller LD. Prognostic and predictive immune gene signatures in breast cancer. Curr Opin Oncol. 2015;27:433–44.View ArticlePubMedGoogle Scholar
- Stefansson OA, Moran S, Gomez A, Sayols S, Arribas-Jorba C, Sandoval J, et al. A DNA methylation-based definition of biologically distinct breast cancer subtypes. Mol Oncol. 2015;9:555–68.View ArticlePubMedGoogle Scholar
- Kamalakaran S, Varadan V, Giercksky Russnes HE, Levy D, Kendall J, Janevski A, et al. DNA methylation patterns in luminal breast cancers differ from non-luminal subtypes and can identify relapse risk independent of other clinical variables. Mol Oncol. 2011;5:77–92.View ArticlePubMedGoogle Scholar
- Nejman D, Straussman R, Steinfeld I, Ruvolo M, Roberts D, Yakhini Z, et al. Molecular rules governing de novo methylation in cancer. Cancer Res. 2014;74:1475–83.View ArticlePubMedGoogle Scholar
- Michaut M, Chin S-F, Majewski I, Severson TM, Bismeijer T, de Koning L, et al. Integration of genomic, transcriptomic and proteomic data identifies two biologically distinct subtypes of invasive lobular breast cancer. Sci Rep. 2016;6:18517.View ArticlePubMedPubMed CentralGoogle Scholar
- Ciriello G, Gatza MLL, Beck AHH, Wilkerson MDD, Rhie SKK, Pastore A, et al. Comprehensive molecular portraits of invasive lobular breast cancer. Cell. 2015;163:506–19.View ArticlePubMedPubMed CentralGoogle Scholar