Abstract
Purpose: High-grade serous ovarian cancers are heterogeneous not only in terms of clinical outcome but also at the molecular level. Our aim was to establish a novel risk classification system based on a gene expression signature for predicting overall survival, leading to suggesting novel therapeutic strategies for high-risk patients.
Experimental Design: In this large-scale cross-platform study of six microarray data sets consisting of 1,054 ovarian cancer patients, we developed a gene expression signature for predicting overall survival by applying elastic net and 10-fold cross-validation to a Japanese data set A (n = 260) and evaluated the signature in five other data sets. Subsequently, we investigated differences in the biological characteristics between high- and low-risk ovarian cancer groups.
Results: An elastic net analysis identified a 126-gene expression signature for predicting overall survival in patients with ovarian cancer using the Japanese data set A (multivariate analysis, P = 4 × 10−20). We validated its predictive ability with five other data sets using multivariate analysis (Tothill's data set, P = 1 × 10−5; Bonome's data set, P = 0.0033; Dressman's data set, P = 0.0016; TCGA data set, P = 0.0027; Japanese data set B, P = 0.021). Through gene ontology and pathway analyses, we identified a significant reduction in expression of immune-response–related genes, especially on the antigen presentation pathway, in high-risk ovarian cancer patients.
Conclusions: This risk classification based on the 126-gene expression signature is an accurate predictor of clinical outcome in patients with advanced stage high-grade serous ovarian cancer and has the potential to develop new therapeutic strategies for high-grade serous ovarian cancer patients. Clin Cancer Res; 18(5); 1374–85. ©2012 AACR.
Using large-scale microarray expression data sets (n = 1,054) by applying an elastic net method, a novel risk classification system for predicting overall survival of patients with advanced stage high-grade serous ovarian cancer based on a 126-gene expression signature was developed and successfully validated. This study has profound significance in clarifying the downregulation of human leukocyte antigen class I antigen presentation machinery that characterizes high-risk ovarian cancer. These results from comprehensive gene expression analysis using large-scale microarray data suggest that our risk classification system might have the potential to optimize treatment of high-grade serous ovarian cancer patients.
Introduction
High-grade serous ovarian cancer comprises approximately 40% of epithelial ovarian cancer cases and is the most aggressive histologic type (1–4). This type of cancer usually presents as advanced stage disease at the time of diagnosis because there are no symptoms present at the early stage and no reliable screening test for early detection (1–4). Patients with advanced stage high-grade serous ovarian cancer generally undergo primary debulking surgery followed by platinum–taxane chemotherapy. However, 30% to 40% of patients recur within 12 months after the standard treatment, and the overall 5-year survival rate remains at approximately 30% (5, 6). Clinicopathologic characteristics, such as the International Federation of Gynecology and Obstetrics (FIGO) stage, histologic grade, and debulking status after primary surgery, are clinically considered important clinical prognostic indicators of ovarian cancer but are insufficient for predicting survival time.
The development of microarray technology has provided new insights into cancer diagnosis and treatment. Large-scale microarray studies in breast cancer have succeeded in clarifying 5 molecular subtypes based on gene expression profiles and in developing genomic biomarker for predicting recurrence in early breast cancer (MammaPrint; refs. 7, 8). Thus, breast cancer treatment strategies are being stratified according to molecular characteristics. In contrast, there are no gene expression signatures with high accuracy and reproducibility for clinical diagnosis and management in patients with ovarian cancer because there is a paucity of ovarian cancer samples available for microarray analysis compared with breast cancer. Although TP53 somatic mutation is present in almost all high-grade serous ovarian cancer and plays an important role in the pathogenesis (9, 10), high-grade serous ovarian cancer exhibits much biological and molecular heterogeneity that should be considered when developing a novel therapeutic strategy for ovarian cancer (10, 11).
Materials and Methods
Clinical samples
Three hundred Japanese patients who were diagnosed with advanced stage high-grade serous ovarian cancer between July 1997 and June 2010 were included in this study. All patients provided written informed consent for the collection of samples and subsequent analysis. Fresh-frozen samples were obtained from primary tumor tissues during debulking surgery prior to chemotherapy. All patients with advanced stage high-grade serous ovarian cancer were treated with platinum–taxane standard chemotherapy after surgery. In principle, patients were seen every 1 to 3 months for the first 2 years. Thereafter, follow-up visits had an interval of 3 to 6 months in the third to fifth year, and 6 to 12 months in the sixth to tenth year. At every follow-up visit, general physical and gynecologic examinations were carried out. CA125 serum levels were routinely determined.
Staging of the disease was assessed according to the criteria of the FIGO (15). Optimal debulking surgery was defined as less than 1 cm of gross residual disease, and suboptimal debulking surgery was defined as more than 1 cm of residual disease. Progression-free survival time was calculated as the interval from primary surgery to disease progression or recurrence. Based on the Response Evaluation Criteria in Solid Tumors (RECIST, version 1.1; 16), disease progression was defined as at least a 20% increase in the sum of the diameters of target lesions, as unequivocal progression of existing nontarget lesions, or as the appearance of one or more new lesions. Overall survival time was calculated as the interval from primary surgery to the death due to ovarian cancer.
The histologic characteristics of surgically resected specimens, which were diagnosed as serous histologic type by 2 pathologists at each institute, were at first assessed on formalin-fixed and paraffin-embedded hematoxylin and eosin sections by 8 gynecologic pathologists who belong to Department of Pathology in 4 institutes (Niigata University, National Defense Medical Collage, Tottori University, Kagoshima City Hospital). Next, these sections were reviewed under central pathologic review by 2 independent gynecologic pathologists (H.T. and T.M.) with no knowledge of patients' clinical data. Histologic subtype was diagnosed based on WHO classification of ovarian tumors (17). The degree of histologic differentiation is determined according to Silverberg classification (18).
Microarray experiments
Frozen tissues containing more than 70% tumor cells upon histologic evaluation were used for RNA extraction. Total RNA was extracted from tissue samples as previously described (19). Five hundred nanograms of total RNA was converted into labeled cRNA with nucleotides coupled to a cyanine 3-CTP (Cy3; PerkinElmer) using the Quick Amp Labeling Kit, one-color (Agilent Technologies). Cy3-labeled cRNA (1.65 μg) was hybridized for 17 hours at 65°C to an Agilent Whole Human Genome Oligo Microarray (G4112F), which carries 60-mer probes to more than 40,000 human transcripts. The hybridized microarray was washed and then scanned in Cy3 channel with the Agilent DNA Microarray Scanner [model G2565CA, n = 260 (Japanese data set A); G2565AA, n = 40 (Japanese data set B)]. Signal intensity per spot was generated from the scanned image with Feature Extraction Software (version 10.1, Japanese data set A; version 9.1, Japanese data set B; Agilent Technologies) with default settings. Spots that did not pass quality control procedures were flagged as “Not Detected.”
Next, we obtained Affymetrix HG-U133Plus2.0 microarray (Affymetrix) data from 10 ovarian cancer samples that had been already been analyzed by the Agilent Whole Human Genome Oligo Microarray. Ten ovarian cancer samples were randomly selected from 260 samples in the Japanese data set A. Microarray experiments were carried out according to the Affymetrix-recommended protocols. Briefly, biotinylated cRNAs were synthesized by GeneChip 3′IVT Express Kit (Affymetrix) from 250 ng total RNA according to the manufacturer's instructions. Biotinylated cRNA yield were checked with the NanoDrop ND-1000 Spectrophotometer. Following fragmentation, 10 μg of cRNA were hybridized for 16 hours at 45°C on GeneChip Human Genome U133 Plus 2.0 Array. GeneChips were washed and stained in the Affymetrix Fluidics Station 450, and scanned with GeneChip Scanner 3000 7G.
The MIAME-compliant microarray data were deposited into the Gene Expression Omnibus data repository (accession number GSE32062 and GSE32063).
Microarray data analysis
We prepared 2 our microarray data sets [Japanese data set A (n = 260) and B (n = 40)] and 4 publicly available large sample-sized (n > 100) microarray data sets [TCGA data set (http://tcga-data.nci.nih.gov/tcga/tcgaHome2.jsp; ref. 11), Tothill's data set (GSE9891; ref. 12), Bonome's data set (GSE26712; ref. 13), and Dressman's data set (14)] to discover predictive biomarkers. Clinical information of publicly available microarray data sets was obtained from their articles and websites. From Tothill's original data set (n = 285), we selected 131 samples that (i) were diagnosed as advanced stage serous adenocarcinoma, (ii) were treated by platinum/taxane–based chemotherapy, and (iii) have clinical data about onset age, stage, grade, surgery, and survival time. Publicly available clinical information in TCGA was downloaded from TCGA Data portal (http://tcga-data.nci.nih.gov/tcga/tcgaHome2.jsp) at June 25, 2011. By the same methods above, 319 samples were selected from 562 microarray data as TCGA data set (Affymetrix HT-HG-U133A). In this study, patients who were treated by adjuvant chemotherapy including molecular-targeted agents were excluded.
On the Agilent platform, data normalization was carried out in GeneSpring GX 11.5 (Agilent Technologies) as follows: (i) threshold raw signals were set to 1.0 and (ii) 75th percentile normalization. Affymetrix microarray data were normalized and summarized with robust multiple average in GeneSpring GX 11.5 (Agilent Technologies). To compare the microarray data sets measured with 4 different platforms (Agilent Whole Human Genome Oligo Microarray, Affymetrix HG-U133A, HG-U133Plus2.0, and HT-HG-U133A), we selected genes common to all platforms based on the Entrez Gene ID and used the Median Rank Score method (20) for cross-platform normalization (Supplementary Fig. S1). Of 22,277 probes that were common among 3 Affymetrix platforms, 20,331 probes were selected to be in one to one relation between probe and gene. Using translation function based on Entrez Gene ID in GeneSpring GX 11.5, the 19,704 transcripts that matched to the 20,331 probes of the Affymetrix platform were extracted from all transcripts on the Agilent platform. Considering differences in microarray platforms, coefficient of correlation (r) in each gene between 2 microarray platforms were measured. Using 10 ovarian cancer microarray data obtained from both Agilent Whole Human Genome Oligo Microarray and Affymetrix HG-U133Plus2.0 (GSE32062), 9,141 genes with high correlation (r > 0.8) were extracted. After Median Rank Score analysis (20), we evaluated median value of each gene in both platforms and selected 3,553 genes with high correlation (r > 0.8) in which absolute value of subtracting median values between 2 platforms was less than 1. Similar analyses were conducted by 16 breast cancer microarray data (GSE17700; ref. 21) from both Affymetrix HG-U133Plus2.0 and HG-U133A. From the resulting 1,746 genes, we removed 60 that were not flagged as “Detected” in more than 90% of the Japanese data set A samples (n = 260), considering them to have either missing or uncertain expression signals. In addition, the data were normalized per gene in each data set by transforming the expression of each gene to obtain a mean of 0 and SD of 1 (Z-transformation) for the cross-platform study.
We analyzed Japanese data set A as a “training set,” Tothill's data set as a “test set,” and the other 4 data sets as “validation sets.” We applied elastic net analysis (22) with the R (23) package glmnet to identify survival-related genes for prediction of prognosis in patients with advanced stage high-grade serous ovarian cancer. Using 10-fold cross-validation, we obtained regression coefficients with optimal penalty parameter for the penalized Cox model, and calculated a prognostic index for each patient as defined by
where βi is the estimated regression coefficient of each gene in the Japanese data set A under elastic net (α = 0.05) and Xi is the Z-transformed expression value of each gene. The estimated regression coefficient of each survival-related gene given by elastic net (α = 0.05) in the Japanese data set A was also applied to calculate a prognostic index for each patient in 5 other data sets using the equation above. We classified all patients into the 2 groups (high- and low-risk groups) by the optimal cutoff value of the prognostic index in the Japanese data set A. Patients were assigned to the “high-risk” group if their prognostic index was more than or equal to cutoff value of prognostic index, whereas “low-risk” group was composed of cases with the prognostic indices that were less than cutoff value. Because risk classification divided by cutoff value 0.1517 indicated a minimum P value (P = 1 × 10−30) when log-rank test was used to compare differences in overall survival between high- and low-risk groups in the Japanese data set A, this value was determined as optimal cutoff value.
Both hierarchical clustering and non-negative factorization (NMF) algorithm (24) were used to assess the similarity of gene expression profiles among 126 survival-related genes. In the 2 major data sets (Japanese data set A and TCGA data set), we constructed a heat map with hierarchical clustering for both samples and genes using Cluster 3.0 (http://bonsai.hgc.jp/~mdehoon/software/cluster/software.htm). Correlation (centered) and complete linkage were selected on similarity metrics and clustering method, respectively. A heat map was visualized with Java TreeView (http://jtreeview.sourceforge.net/). NMF-consensus matrices averaging 50 connectivity matrices were computed at K = 2–7 (as the number of subclasses modeled) for 126 genes. With genes appearing along both the horizontal and vertical axes of the consensus matrices, consistency in the gene-pair clusters was visualized. We determined the optimal number of gene clusters by the cophenetic correlation coefficient that provides a scalar summary of global clustering robustness across the consensus matrix in the Japanese data set A and TCGA data set.
We conducted a volcano plot analysis to extract differentially expressed transcripts between high- and low-risk ovarian cancer groups with the 2 major data sets [Japanese data set A (n = 260) and TCGA data set (n = 319)]. When the volcano plot analysis was conducted in GeneSpring GX 11.5 (Agilent Technologies), we used gene expression data prior to Z-transformation normalization (Fig. 3).
To investigate the biological functions of gene expression signatures, we used GO Ontology Browser, embedded in GeneSpring GX11.5 (Agilent Technologies). The GO Ontology Browser was used to analyze which categories of gene ontology were statistically overrepresented among the gene list obtained. Statistical significance was determined by Fisher exact test, followed by multiple testing corrections by the Benjamini and Yekutieli false discovery rate method (25) To avoid bias of gene extraction by volcano plot analysis in GO analysis, Gene Set Enrichment Analysis (GSEA; ref. 26) was conducted with genes prior to gene selection by volcano plot analysis. Analysis settings in GSEA (software version 2.07) were as follows: (i) gene sets database: c5.all.v2.5.symbols.gmt [gene ontology], (ii) number of permutations: 1,000, (iii) collapse data set to gene symbol: true, (iv) permutation type: phenotype, (v) chip platform(s): Agilent_HumanGenome.chip or HT_HG_U133A.chip, and (vi) other settings: default.
Furthermore, we used the Core Analysis tool in the Ingenuity Pathway Analysis (IPA) system to analyze networks and pathways for a set of genes. Q value was calculated by Fisher exact test with the Benjamini and Hochberg correction (27). Q < 0.25 was considered as significant in GO, GSEA, and IPA.
Single-nucleotide polymorphism array experiments
We isolated genomic DNA from tumor tissues with a phenol-chloroform extraction method and from normal lymphocytes using the QIAamp DNA Blood Maxi Kit (QIAGEN). Single-nucleotide polymorphism (SNP) array experiments with the Genome-Wide Human SNP Array 6.0 (Affymetrix) were carried out at Niigata University (details in the Supplementary Methods). SNP array data were analyzed by the Partek Genomic Suite 6.5 (Partek Inc.) to investigate copy number variations in 30 genes involved in the antigen presentation pathway.
Immunoshitochemical analysis
A monoclonal anti-human CD8 antibody (M7103; 1:100; DAKO) was used for immunohistochemical staining of CD8 T lymphocyte in tumor tissues. Details of the immunohistochemistry method and sample selections are described in the Supplementary Methods.
Statistical analysis
Standard statistical tests including Pearson correlation analysis, unpaired t tests, 1-way ANOVA, Fisher exact tests, log-rank tests, and Cox proportional hazard model analysis were used to analyze the clinical data, as appropriate. Analyses of clinical data were conducted with JMP (version 8; SAS Institute) and GraphPad PRISM (version 4.0; GraphPad Software).
Results
The distribution of the clinical variables for each microarray data set is shown in Table 1. To compare the microarray data sets measured with 4 different platforms, 1,686 genes were selected (Supplementary Fig. S1). An elastic net analysis (22) using 10-fold cross-validation on Japanese data set A (training set, n = 260) identified a 126-gene signature for predicting overall survival in patents with advanced stage high-grade serous ovarian cancer (Supplementary Table S1). After calculating the prognostic index for each sample from the 126-gene expression signature as reported previously (19), we divided the training set into high- and low-risk groups based on the optimal cutoff value (0.1517) of the prognostic index and verified the high predictive power of this risk classification (log-rank P = 1 × 10−30; Multivariate Cox P = 4 × 10−20; HR = 6.203; 95% CI = 4.239–9.123; Fig. 1A and Table 2).
Data set . | Japanese data set A . | Tothill's data seta . | Bonome's data seta . | Dressman's data seta . | TCGA data setc . | Japanese data set B . |
---|---|---|---|---|---|---|
Number | 260 | 131 | 185 | 119 | 319 | 40 |
Age | 58.2 ± 10.8 | 58.4 ± 9.8 | 62 ± 12 | N/Ab | 59.5 ± 11.3 | 56.2 ± 9.6 |
Histology | ||||||
Serous | 260 | 131 | 166 | 119 | 319 | 40 |
Others | 0 | 0 | 19 | 0 | 0 | 0 |
Stage | ||||||
III | 204 | 123 | 144 | 99 | 267 | 31 |
IV | 56 | 8 | 41 | 20 | 52 | 9 |
Grade | Silverberg | Silverberg | N/Ab | N/Ab | N/Ab | Silverberg |
1 | 0 | 0 | 0 | 3 | 0 | 0 |
2 | 131 | 51 | 40 | 57 | 29 | 23 |
3 | 129 | 80 | 144 | 59 | 289 | 17 |
4 | — | — | 3 | — | 1 | — |
Surgery status | ||||||
Optimal | 103 | 78 | 92 | 63 | 236 | 19 |
Suboptimal | 157 | 45 | 93 | 56 | 83 | 21 |
Chemotherapy | ||||||
Platinum | 260 | 131 | 185 | 119 | 319 | 40 |
Taxane | 260 | 131 | N/A | 82 | 319 | 40 |
Follow-up period (mo) | ||||||
Median | 42 | 29 | 38 | 34 | 31 | 39 |
Range | 1–128 | 6–79 | 1–164 | 1–185 | 1–154 | 7–111 |
Number of deaths | 121 | 60 | 129 | 69 | 187 | 22 |
Median survival (mo) | 60 | 44 | 46 | 69 | 42 | 54 |
Microarray platform | Agilent Whole Human Genome Oligo Microarray | Affymetrix HG-U133 Plus2.0 | Affymetrix HG-U133A | Affymetrix HG-U133A | Affymetix HT-HG-U133A | Agilent Whole Human Genome Oligo Microarray |
Data set . | Japanese data set A . | Tothill's data seta . | Bonome's data seta . | Dressman's data seta . | TCGA data setc . | Japanese data set B . |
---|---|---|---|---|---|---|
Number | 260 | 131 | 185 | 119 | 319 | 40 |
Age | 58.2 ± 10.8 | 58.4 ± 9.8 | 62 ± 12 | N/Ab | 59.5 ± 11.3 | 56.2 ± 9.6 |
Histology | ||||||
Serous | 260 | 131 | 166 | 119 | 319 | 40 |
Others | 0 | 0 | 19 | 0 | 0 | 0 |
Stage | ||||||
III | 204 | 123 | 144 | 99 | 267 | 31 |
IV | 56 | 8 | 41 | 20 | 52 | 9 |
Grade | Silverberg | Silverberg | N/Ab | N/Ab | N/Ab | Silverberg |
1 | 0 | 0 | 0 | 3 | 0 | 0 |
2 | 131 | 51 | 40 | 57 | 29 | 23 |
3 | 129 | 80 | 144 | 59 | 289 | 17 |
4 | — | — | 3 | — | 1 | — |
Surgery status | ||||||
Optimal | 103 | 78 | 92 | 63 | 236 | 19 |
Suboptimal | 157 | 45 | 93 | 56 | 83 | 21 |
Chemotherapy | ||||||
Platinum | 260 | 131 | 185 | 119 | 319 | 40 |
Taxane | 260 | 131 | N/A | 82 | 319 | 40 |
Follow-up period (mo) | ||||||
Median | 42 | 29 | 38 | 34 | 31 | 39 |
Range | 1–128 | 6–79 | 1–164 | 1–185 | 1–154 | 7–111 |
Number of deaths | 121 | 60 | 129 | 69 | 187 | 22 |
Median survival (mo) | 60 | 44 | 46 | 69 | 42 | 54 |
Microarray platform | Agilent Whole Human Genome Oligo Microarray | Affymetrix HG-U133 Plus2.0 | Affymetrix HG-U133A | Affymetrix HG-U133A | Affymetix HT-HG-U133A | Agilent Whole Human Genome Oligo Microarray |
aClinical information in these 3 data sets were obtained from their articles and website.
bNo available information was described as N/A.
cPublicly available clinical information in TCGA was downloaded from TCGA Data portal (http://tcga-data.nci.nih.gov/tcga/tcgaHome2.jsp) at June 25, 2011.
. | Univariate analysis . | P . | Multivariate analysis . | P . |
---|---|---|---|---|
. | HR (95% CI) . | . | HR (95% CI) . | . |
Japanese data set A | ||||
Age | 1.011 (0.993–1.029) | 0.21 | 1.006 (0.988–1.024) | 0.52 |
Stage IV (vs. stage III) | 1.465 (0.968–2.165) | 0.07 | 1.340 (0.884–1.984) | 0.16 |
Optimal surgery (vs. suboptimal) | 0.499 (0.334–0.730) | 0.0003 | 0.643 (0.428–0.949) | 0.026 |
High (vs. low) | 6.823 (4.692–9.972) | 2 × 10−22 | 6.203 (4.239–9.123) | 4 × 10−20 |
Tothill's data set | ||||
Age | 1.005 (0.977–1.035) | 0.73 | 1.006 (0.978–1.035) | 0.67 |
Stage IV (vs. stage III) | 2.264 (0.785–5.176) | 0.12 | 2.951 (0.983–7.235) | 0.053 |
Optimal surgery (vs. suboptimal) | 0.887 (0.525–1.517) | 0.66 | 0.910 (0.528–1.592) | 0.74 |
High (vs. low) | 3.443 (1.967–5.962) | 3 × 10−5 | 3.728 (2.110–6.532) | 1 × 10−5 |
Bonome's data set | ||||
Optimal surgery (vs. suboptimal) | 0.592 (0.414–0.840) | 0.0032 | 0.628 (0.438–0.894) | 0.0097 |
High (vs. low) | 2.034 (1.341–3.012) | 0.0011 | 1.897 (1.247–2.818) | 0.0033 |
Dressman's data set | ||||
Optimal surgery (vs. suboptimal) | 0.685 (0.425–1.102) | 0.12 | 0.602 (0.370–0.977) | 0.04 |
High (vs. low) | 2.397 (1.335–4.140) | 0.0042 | 2.687 (1.480–4.705) | 0.0016 |
TCGA data set | ||||
Age | 1.014 (1.001–1.028) | 0.038 | 1.011 (0.997–1.025) | 0.13 |
Stage IV (vs. stage III) | 1.018 (0.683–1.469) | 0.93 | 1.034 (0.692–1.499) | 0.86 |
Optimal surgery (vs. suboptimal) | 0.909 (0.668–1.254) | 0.56 | 0.903 (0.654–1.261) | 0.54 |
High (vs. low) | 1.712 (1.234–2.345) | 0.0015 | 1.680 (1.202–2.321) | 0.0027 |
Japanese data set B | ||||
Age | 1.059 (1.007–1.118) | 0.024 | 1.058 (0.999–1.131) | 0.057 |
Stage IV (vs. stage III) | 1.696 (0.542–4.501) | 0.34 | 2.998 (0.808–10.60) | 0.098 |
Optimal surgery (vs. suboptimal) | 0.582 (0.231–1.364) | 0.22 | 0.510 (0.175–1.413) | 0.20 |
High (vs. low) | 4.100 (1.376–11.25) | 0.013 | 4.468 (1.265–15.49) | 0.021 |
. | Univariate analysis . | P . | Multivariate analysis . | P . |
---|---|---|---|---|
. | HR (95% CI) . | . | HR (95% CI) . | . |
Japanese data set A | ||||
Age | 1.011 (0.993–1.029) | 0.21 | 1.006 (0.988–1.024) | 0.52 |
Stage IV (vs. stage III) | 1.465 (0.968–2.165) | 0.07 | 1.340 (0.884–1.984) | 0.16 |
Optimal surgery (vs. suboptimal) | 0.499 (0.334–0.730) | 0.0003 | 0.643 (0.428–0.949) | 0.026 |
High (vs. low) | 6.823 (4.692–9.972) | 2 × 10−22 | 6.203 (4.239–9.123) | 4 × 10−20 |
Tothill's data set | ||||
Age | 1.005 (0.977–1.035) | 0.73 | 1.006 (0.978–1.035) | 0.67 |
Stage IV (vs. stage III) | 2.264 (0.785–5.176) | 0.12 | 2.951 (0.983–7.235) | 0.053 |
Optimal surgery (vs. suboptimal) | 0.887 (0.525–1.517) | 0.66 | 0.910 (0.528–1.592) | 0.74 |
High (vs. low) | 3.443 (1.967–5.962) | 3 × 10−5 | 3.728 (2.110–6.532) | 1 × 10−5 |
Bonome's data set | ||||
Optimal surgery (vs. suboptimal) | 0.592 (0.414–0.840) | 0.0032 | 0.628 (0.438–0.894) | 0.0097 |
High (vs. low) | 2.034 (1.341–3.012) | 0.0011 | 1.897 (1.247–2.818) | 0.0033 |
Dressman's data set | ||||
Optimal surgery (vs. suboptimal) | 0.685 (0.425–1.102) | 0.12 | 0.602 (0.370–0.977) | 0.04 |
High (vs. low) | 2.397 (1.335–4.140) | 0.0042 | 2.687 (1.480–4.705) | 0.0016 |
TCGA data set | ||||
Age | 1.014 (1.001–1.028) | 0.038 | 1.011 (0.997–1.025) | 0.13 |
Stage IV (vs. stage III) | 1.018 (0.683–1.469) | 0.93 | 1.034 (0.692–1.499) | 0.86 |
Optimal surgery (vs. suboptimal) | 0.909 (0.668–1.254) | 0.56 | 0.903 (0.654–1.261) | 0.54 |
High (vs. low) | 1.712 (1.234–2.345) | 0.0015 | 1.680 (1.202–2.321) | 0.0027 |
Japanese data set B | ||||
Age | 1.059 (1.007–1.118) | 0.024 | 1.058 (0.999–1.131) | 0.057 |
Stage IV (vs. stage III) | 1.696 (0.542–4.501) | 0.34 | 2.998 (0.808–10.60) | 0.098 |
Optimal surgery (vs. suboptimal) | 0.582 (0.231–1.364) | 0.22 | 0.510 (0.175–1.413) | 0.20 |
High (vs. low) | 4.100 (1.376–11.25) | 0.013 | 4.468 (1.265–15.49) | 0.021 |
The predictive power of the 126-gene signature was tested with Tothill's data set (n = 131; ref. 12). Both Kaplan–Meier survival and multivariate analysis showed that this risk classification was significantly associated with overall survival time in Tothill's data set (log-rank P = 3 × 10−6; Multivariate Cox P = 1 × 10−5; HR = 3.728; 95% CI = 2.110–6.532; Fig. 1B and Table 2). Next, we assessed the predictive power of the 126-gene expression signature on the 3 data sets in which microarray data were obtained from Affymetrix platforms and confirmed that this risk classification indicated similar results in the 3 data sets (Fig. 1C–E and Table 2). Furthermore, to exclude the influence of differences in microarray platforms, we prepared Japanese data set B from the same platform as the training set. Despite the small sample size, our risk classification showed a significant association with overall survival time in the Japanese data set B (Fig. 1F and Table 2). This risk classification showed high predictive accuracy even though overall survival was capped at 60 months (Supplementary Fig. S2; ref. 10). Moreover, the high-risk groups had shorter progression-free survival times compared with the low-risk groups in 4 microarray data sets with available progression-free survival times (Supplementary Fig. S3). In addition, we applied 193 overall survival related gene signature that have been recently reported by TCGA research networks (10) to our large data set (Japanese data set A), and 191 of 193 genes were available in our data set. High-risk ovarian cancer patients based on 191-gene signature model had significantly poor prognosis compared with low-risk patients (Supplementary Methods and Supplementary Fig. S4).
We next investigated the molecular characteristics of the 126-gene signature with both hierarchical clustering and NMF (24) analyses using the 2 major data sets (Japanese data set A and TCGA data set). Hierarchical clustering of the expression of 126 genes in the Japanese data set A (n = 260) resulted in 2 gene clusters (cluster 1 and 2; Fig. 2A). Cluster 1 (outlined in yellow in Fig. 2A) comprised 23 genes, 22 of which were common genes in a small cluster of NMF class assignment for K = 2 (Supplementary Fig. S5). Similarly, this gene cluster was seen both in hierarchical clustering and NMF analyses of TCGA data set, and the 20 genes in TCGA's cluster were consistent with these 23 genes in cluster 1 (Fig. 2B and Supplementary Fig. S5). GO analysis of these 23 genes revealed 132 significantly overrepresented GO categories; the top 20 GO categories are shown in the Supplementary Table S2. Specifically, immunity-related categories were enriched in these 23 genes and were downregulated in the high-risk groups compared with the low-risk groups. On the contrary, only 1 category [cytoplasm (GO0005737)] belonging to “cellular component” was significantly overrepresented in 103 genes of cluster 2.
To clarify the biological significance of our risk classification based on the 126-gene signature in advanced stage high-grade serous ovarian cancer, we examined differences in molecular biological characteristics between high- and low-risk groups. For subsequent analyses, all transcripts on each platform were used to avoid the previous gene number limitations, which were for the cross-platform study. We extracted the genes differentially expressed between the 2 groups by carrying out a volcano plot analysis with the 2 major data sets (Japanese data set A and TCGA data set; Supplementary Table S3). Figure 3A shows that 1,109 and 1,381 transcripts were differentially expressed between the high- and low-risk groups in the Japanese data set A and TCGA data set, respectively. GO analysis of these transcripts indicated that 9 out of 20 top-ranked GO categories were common between both data sets and involved in the immune system (Table 3 and 4).
. | Genes within GO category . | . | |
---|---|---|---|
GO categorya . | Number . | Percentage . | −Log10Qb . |
Immune system process (GO:0002376) | 206 | 29.0 | 45.0 |
Immune response (GO:0006955) | 172 | 24.2 | 45.0 |
Defense response (GO:0006952, 0002217, 0042829) | 110 | 15.5 | 40.3 |
Response to stimulus (GO:0050896, 0051869) | 251 | 35.4 | 29.2 |
Inflammatory response (GO:0006954) | 68 | 9.6 | 26.2 |
Positive regulation of immune system process (GO:0002684) | 40 | 5.6 | 25.2 |
Regulation of immune system process (GO:0002682) | 52 | 7.3 | 24.5 |
Antigen processing and presentation (GO:0019882,0030333) | 31 | 4.4 | 22.7 |
Signal transducer activity (GO:0004871, 0005062, 0009369, 0009370) | 169 | 23.8 | 22.1 |
Molecular transducer activity (GO:0060089) | 169 | 23.8 | 22.1 |
Signal transduction (GO:0007165) | 224 | 31.5 | 20.8 |
Leukocyte activation (GO:0045321) | 47 | 6.6 | 20.1 |
Cell activation (GO:0001775) | 48 | 6.8 | 19.7 |
Regulation of immune response (GO:0050776) | 32 | 4.5 | 18.9 |
Response to wounding (GO:0009611, 0002245) | 70 | 9.9 | 18.7 |
Signal transmission (GO:0023060) | 224 | 31.5 | 18.2 |
Signaling process (GO:0023046) | 224 | 31.5 | 18.2 |
Regulation of cell activation (GO:0050865) | 30 | 4.2 | 17.2 |
Positive regulation of cell activation (GO:0050867) | 24 | 3.4 | 17.0 |
Regulation of lymphocyte activation (GO:0051249) | 29 | 4.1 | 17.0 |
. | Genes within GO category . | . | |
---|---|---|---|
GO categorya . | Number . | Percentage . | −Log10Qb . |
Immune system process (GO:0002376) | 206 | 29.0 | 45.0 |
Immune response (GO:0006955) | 172 | 24.2 | 45.0 |
Defense response (GO:0006952, 0002217, 0042829) | 110 | 15.5 | 40.3 |
Response to stimulus (GO:0050896, 0051869) | 251 | 35.4 | 29.2 |
Inflammatory response (GO:0006954) | 68 | 9.6 | 26.2 |
Positive regulation of immune system process (GO:0002684) | 40 | 5.6 | 25.2 |
Regulation of immune system process (GO:0002682) | 52 | 7.3 | 24.5 |
Antigen processing and presentation (GO:0019882,0030333) | 31 | 4.4 | 22.7 |
Signal transducer activity (GO:0004871, 0005062, 0009369, 0009370) | 169 | 23.8 | 22.1 |
Molecular transducer activity (GO:0060089) | 169 | 23.8 | 22.1 |
Signal transduction (GO:0007165) | 224 | 31.5 | 20.8 |
Leukocyte activation (GO:0045321) | 47 | 6.6 | 20.1 |
Cell activation (GO:0001775) | 48 | 6.8 | 19.7 |
Regulation of immune response (GO:0050776) | 32 | 4.5 | 18.9 |
Response to wounding (GO:0009611, 0002245) | 70 | 9.9 | 18.7 |
Signal transmission (GO:0023060) | 224 | 31.5 | 18.2 |
Signaling process (GO:0023046) | 224 | 31.5 | 18.2 |
Regulation of cell activation (GO:0050865) | 30 | 4.2 | 17.2 |
Positive regulation of cell activation (GO:0050867) | 24 | 3.4 | 17.0 |
Regulation of lymphocyte activation (GO:0051249) | 29 | 4.1 | 17.0 |
aBold font denotes common categories included in top 20 lists both Japanese data set A and TCGA data set.
bQ value was determined by Fisher's exact test with Benjamini–Yekutieli correction.
. | Genes within GO category . | . | |
---|---|---|---|
GO categorya . | Number . | Percentage . | −Log10Qb . |
Immune system process (GO:0002376) | 129 | 20.5 | 20.9 |
Immune response (GO:0006955) | 115 | 18.3 | 20.8 |
Defense response (GO:0006952, 0002217, 0042829) | 78 | 12.4 | 10.3 |
Antigen processing and presentation (GO:0019882, 0030333) | 25 | 4.0 | 8.4 |
Inflammatory response (GO:0006954) | 50 | 7.9 | 8.2 |
Antigen processing and presentation of peptide antigen (GO:0048002) | 13 | 2.1 | 6.9 |
Antigen processing and presentation of exogenous peptide antigen (GO:0002478) | 6 | 1.0 | 5.2 |
MHC class I peptide loading complex (GO:0042824) | 9 | 1.4 | 5.2 |
MHC protein complex (GO:0042611) | 18 | 2.9 | 5.2 |
TAP complex (GO:00042825) | 8 | 1.3 | 5.2 |
MHC protein binding (GO:0042287) | 12 | 1.9 | 5.2 |
Antigen processing and presentation of exogenous antigen (GO:0019884) | 7 | 1.1 | 5.0 |
Response to stimulus (GO:0050896, 0051869) | 183 | 29.1 | 4.4 |
Antigen processing and presentation of peptide antigen via MHC class I (GO:0002474) | 7 | 1.1 | 3.9 |
Antigen processing and presentation of peptide or polysaccharide antigen via MHC class II (GO:0002504) | 12 | 1.9 | 3.9 |
Regulation of immune response (GO:0050776) | 6 | 1.0 | 3.8 |
MHC class I protein binding (GO:0042288) | 10 | 1.6 | 3.8 |
Response to wounding (GO:0009611, 0002245) | 50 | 7.9 | 3.8 |
Positive regulation of immune response (GO:0050778) | 6 | 1.0 | 3.8 |
Positive regulation of immune system process (GO:0002684) | 6 | 1.0 | 3.8 |
. | Genes within GO category . | . | |
---|---|---|---|
GO categorya . | Number . | Percentage . | −Log10Qb . |
Immune system process (GO:0002376) | 129 | 20.5 | 20.9 |
Immune response (GO:0006955) | 115 | 18.3 | 20.8 |
Defense response (GO:0006952, 0002217, 0042829) | 78 | 12.4 | 10.3 |
Antigen processing and presentation (GO:0019882, 0030333) | 25 | 4.0 | 8.4 |
Inflammatory response (GO:0006954) | 50 | 7.9 | 8.2 |
Antigen processing and presentation of peptide antigen (GO:0048002) | 13 | 2.1 | 6.9 |
Antigen processing and presentation of exogenous peptide antigen (GO:0002478) | 6 | 1.0 | 5.2 |
MHC class I peptide loading complex (GO:0042824) | 9 | 1.4 | 5.2 |
MHC protein complex (GO:0042611) | 18 | 2.9 | 5.2 |
TAP complex (GO:00042825) | 8 | 1.3 | 5.2 |
MHC protein binding (GO:0042287) | 12 | 1.9 | 5.2 |
Antigen processing and presentation of exogenous antigen (GO:0019884) | 7 | 1.1 | 5.0 |
Response to stimulus (GO:0050896, 0051869) | 183 | 29.1 | 4.4 |
Antigen processing and presentation of peptide antigen via MHC class I (GO:0002474) | 7 | 1.1 | 3.9 |
Antigen processing and presentation of peptide or polysaccharide antigen via MHC class II (GO:0002504) | 12 | 1.9 | 3.9 |
Regulation of immune response (GO:0050776) | 6 | 1.0 | 3.8 |
MHC class I protein binding (GO:0042288) | 10 | 1.6 | 3.8 |
Response to wounding (GO:0009611, 0002245) | 50 | 7.9 | 3.8 |
Positive regulation of immune response (GO:0050778) | 6 | 1.0 | 3.8 |
Positive regulation of immune system process (GO:0002684) | 6 | 1.0 | 3.8 |
aBold font denotes common categories included in top 20 lists both Japanese data set A and TCGA data set.
bQ value was determined by Fisher's exact test with Benjamini–Yekutieli correction.
To exclude the influence of cutoff in fold change or P value in volcano plot analysis, we reevaluated differences in the biological characteristics between high- and low-risk groups using GSEA (26). Four immunity-related GO categories (immune system response, immune response, defense response, and inflammatory response) were also included in the list of 20 top-ranked GO categories when GSEA was carried out with 21,785 (Japanese data set A) and 14,023 (TCGA data set) transcripts prior to gene extraction (Supplementary Table S4).
Furthermore, possible functional relations among differentially expressed genes between high- and low-groups were investigated with pathway analysis. Of 20 top-ranked pathways notably enriched in 1,109 transcripts of the Japanese data set A and 1,381 transcripts of the TCGA data set, 14 pathways were common between both data sets (Supplementary Table S5). In particular, the antigen presentation pathway (Fig. 3B) was the most significantly overrepresented pathway in the 2 data sets (Japanese data set A, Q = 2.0 × 10−29; TCGA data set, Q = 5.0 × 10−14). In the antigen presentation pathway, 30 transcripts in this pathway were significantly downregulated in the high-risk group (Fig. 3B and Supplementary Table S6). We further examined whether molecular defects of human leukocyte antigen (HLA) class I antigen presentation machinery components were caused by structural alterations using SNP array data. Our SNP array data showed that genes of HLA class I antigen presentation machinery were not deleted (Supplementary Table S6). In the TCGA data set, only 10% of cases (24 of 242) had deletions in HLA class I genes. On the other hand, 89% of high-risk cases (55 of 62) without deletion in HLA class I genes showed significantly lower expressions in HLA class I genes compared with those in low-risk groups (Supplementary Fig. S6).
On the basis of this result, we assessed the status of tumor-infiltrating lymphocyte reflecting an immune response against ovarian cancer. By immunohistochemically staining for CD8-positive T lymphocytes in a subset of the Japanese data set (n = 30), we revealed that the number of CD8 T lymphocytes infiltrating into tumor tissues was significantly decreased in the high-risk group (n = 10) compared with the low-risk group (n = 20) and was clearly correlated with the expression level of HLA class I genes (Supplementary Fig. S7).
Discussion
In this study, we established a novel risk classification system based on the 126-gene expression signature for predicting overall survival time in patients with advanced stage high-grade serous ovarian cancer. The significant association between our risk classification and overall survival time was indicated among the 6 microarray data sets.
In expression microarray analyses, there is a well-known “curse of dimensionality” problem that the number of genes is much larger than the number of samples. The curse of dimensionality leads to a concern about the reliability of the selected genes or an overfitting phenomenon, and using a large-scale microarray data is a simple and effective method to overcome this problem. From this theory, we planned a large-scale cross-platform study using 4 publicly available microarray data sets consisting of more than 100 samples. To carry out a cross-platform study, we took considerable care to reduce the influence of differences in microarray platforms with the Median Rank Score method and matching probes based on highly correlated expression (r > 0.8) between the 2 platforms. As a result, this risk classification system showed a high predictive ability in the 4 microarray data sets from different platforms. Our risk classification was suitably fitting to Japanese data set B despite the small sample size. This might reflect the same platform in the microarray experiments and similarity in clinical backgrounds between Japanese data set A and B. Intriguingly, 193 survival-related gene signature that has been reported by TCGA research networks (10) was also significantly associated with overall survival in our large data set (Supplementary Fig. S4), and 5 genes (B4GALT5, CXCL9, ID4, MTRF1, and SLC7A11) were common between their and our gene signatures despite different analytic processes (Supplementary Table 2). These genes might contribute strongly to developing risk classification system for high-grade serous ovarian cancer patients.
We ascertained that alterations to the immune system in cancer cells are one of the most important factors affecting survival of patients with advanced stage, high-grade serous ovarian cancer and that high-risk ovarian cancer was well characterized by the downregulation of the antigen presentation pathway at the molecular level.
The 2 large-scale subclassification studies (10, 11) based on gene expression analyses have indicated that there are at least 4 subclasses (immunoreactive, differentiated, proliferative, and mesenchymal) in high-grade serous ovarian cancer. Although there was no statistically significant difference in clinical outcome among 4 subclasses, this finding that immunity-related genes are identified in large-scale gene expression profiles despite different analytic approaches suggests that alteration of immune activity might play an important role in the molecular characterization of ovarian cancer.
Previous findings (28–30) that the presence of tumor-infiltrating T lymphocytes is associated with long survival in ovarian cancer patients are consistent with the results from our comprehensive gene expression analysis using large-scale microarray data, suggests that the presence of immune responses to cancer cells would influence the biological phenotype of ovarian cancer. Although several antigen-specific active immunotherapy studies have been conducted, the clinical benefit of this approach has not yet been shown in large randomized-controlled trials (31, 32). Previous reports (28, 33, 34) and our data indicated that defects in HLA class I antigen presentation machinery would decrease recruitment of tumor-infiltrating lymphocytes, leading to poor prognosis in cancer patients because of a reduction in antitumor immune activity. Therefore, upregulation of HLA class I antigen presentation pathway might be one of efficient therapeutic approach for patients with high-risk serous ovarian cancer, especially when antigen-specific active immunotherapy is selected as a therapeutic strategy (35).
Defects of HLA class I gene expression occurs at the genetic, epigenetic, transcriptional, and posttranscriptional levels, and are classified into the 2 main groups: reversible and irreversible defects (36). Our data show that the frequency of deletion leading to irreversible defect of HLA class I gene expression was low in high-grade serous ovarian cancer. Only a few HLA class I gene mutations have been described thus far (36). Downregulation of antigen presentation machinery components such as TAP1/2 or B2M also result in reversible defects in HLA class I molecules (37). Interestingly, the expression levels of genes in the antigen presentation pathway were positively correlated in the Japanese data set A (Supplementary Fig. S8). IFN-γ or other cytokines stimulation induces the expression of genes in the antigen presentation pathway (38). Moreover, histone deacetylase (HDAC) inhibitors increase the expression of genes in the antigen presentation pathway such as TAP1/2 (39, 40) and B2M (41) in several cancer cells leading to upregulating the antigen presentation pathway. In addition, HDAC inhibitors that are promising anticancer drugs show a synergistic effect with taxane or platinum drugs, which were used as standard adjuvant chemotherapy for ovarian cancer, both in vitro and in vivo (42, 43). These genetic and epigenetic activations of the antigen presentation pathway might induce immune recognition of ovarian cancer cells and enhance antitumor immune responses.
In summary, our data suggest that this predictive biomarker based on the 126-gene signature could identify patients who should not expect long-term survival by standard treatment and that activation of the antigen presentation pathway in tumor cells is an important key in new therapeutic strategies for ovarian cancer.
Disclosure of Potential Conflicts of Interest
No potential conflicts of interest were disclosed.
Acknowledgments
The authors thank Prof. Y. Nakamura for his help, the tissue donors and supporting medical staff for making this study possible, C. Seki and A. Yukawa for their technical assistance, T. Mizuochi for discussion, Prof. S. G. Silverberg for his enthusiastic help with the pathologic review, and K. Boroevich for helpful comments on the manuscript.
The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.