Abstract
Background: While numerous susceptibility loci for epithelial ovarian cancer (EOC) have been identified, few associations have been reported with overall survival. In the absence of common prognostic genetic markers, we hypothesize that rare coding variants may be associated with overall EOC survival and assessed their contribution in two exome-based genotyping projects of the Ovarian Cancer Association Consortium (OCAC).
Methods: The primary patient set (Set 1) included 14 independent EOC studies (4,293 patients) and 227,892 variants, and a secondary patient set (Set 2) included six additional EOC studies (1,744 patients) and 114,620 variants. Because power to detect rare variants individually is reduced, gene-level tests were conducted. Sets were analyzed separately at individual variants and by gene, and then combined with meta-analyses (73,203 variants and 13,163 genes overlapped).
Results: No individual variant reached genome-wide statistical significance. A SNP previously implicated to be associated with EOC risk and, to a lesser extent, survival, rs8170, showed the strongest evidence of association with survival and similar effect size estimates across sets (Pmeta = 1.1E−6, HRSet1 = 1.17, HRSet2 = 1.14). Rare variants in ATG2B, an autophagy gene important for apoptosis, were significantly associated with survival after multiple testing correction (Pmeta = 1.1E−6; Pcorrected = 0.01).
Conclusions: Common variant rs8170 and rare variants in ATG2B may be associated with EOC overall survival, although further study is needed.
Impact: This study represents the first exome-wide association study of EOC survival to include rare variant analyses, and suggests that complementary single variant and gene-level analyses in large studies are needed to identify rare variants that warrant follow-up study. Cancer Epidemiol Biomarkers Prev; 25(3); 446–54. ©2016 AACR.
This article is featured in Highlights of This Issue, p. 427
Introduction
Epithelial ovarian cancer (EOC) is diagnosed in over 230,000 women worldwide every year and is a leading cause of cancer-related death (1). Most women are diagnosed with advanced stage disease, when five-year survival rates are low (2, 3). Several clinical and demographic factors are associated with survival, such as stage, grade, and histologic subtype (4), but few germline prognostic variants have been identified (5). The strongest known genetic risk factors for EOC are rare BRCA1 and BRCA2 mutations (6), which also confer differences in clinical characteristics including improved five-year, but not ten-year, survival (7–9). Although genome-wide association studies (GWAS) to identify common risk loci have been fruitful (10–15), studies of survival have been less so (5). For example, variant rs8170 at the 19p13 locus, a region associated with ovarian and breast cancer risk phenotypes and known to interact with BRCA1 (10, 16), was associated with EOC survival in an initial phase of a GWAS (HR = 1.35, P = 2.4E−4) conducted by the Ovarian Cancer Association Consortium (OCAC), but not replicated in an independent dataset (HR = 1.01, P = 0.85) (10).
While the lack of identified common germline prognostic markers from earlier GWAS may be due to inadequate power (due to limited sample size or study heterogeneity), it is also possible that rare variants (in addition to those in BRCA1 and BRCA2) not captured in GWAS arrays may be associated with overall EOC survival. In fact, it has been suggested that multiple rare variants of large effect could collectively be responsible for some of the “missing heritability” not explained by the common variants of modest effect identified through GWAS (17). Motivated by this rare variant hypothesis (18), commercial genotype arrays based on exome sequencing studies that attempted to capture all variants (common and rare) within coding regions have emerged as a new approach. As single marker tests have very little power to detect association with rare variants, approaches that pool information across all variants within a gene region may provide improved power (19–23). We hypothesize that rare variants, either individually or collectively across a gene, may be associated with overall survival in EOC. To test this hypothesis, we combined data from two exome-based OCAC genotyping projects that used commercial arrays based on 16 major exome-sequencing projects (24). With over 6,000 EOC patients, this is the first exome-wide rare variant assessment of genetic associations with EOC overall survival.
Materials and Methods
Study participants
Study participants included 20 independent studies of EOC (Supplementary Table S1), which consisted of two sample sets that were genotyped separately on different platforms.
Set 1.
For Set 1, 14 independent studies of EOC (DOV, HAW, HOP, LAX, MAC, MAY, MSK, NCO, NEC, NJO, ORE, POL, UCI, USC) served as the primary sample (Supplementary Table S1), including 6,293 patients. Patients consisted of women aged 18 and older with a pathologically confirmed primary invasive EOC, fallopian tube cancer, or primary peritoneal cancer, excluding Brenner tumors or those missing tumor histology. Patients were excluded with incomplete survival time information, such as missing vital status or time from diagnosis to follow-up, resulting in 4,976 patients. To avoid potential issues with population stratification, only patients of European ancestry were analyzed (N = 4,293; 2,257 deaths at 10 years postdiagnosis).
Set 2.
The smaller Set 2 consisted of six additional independent studies of invasive EOC (AUS, MAL, POC, RMH, SEA, UKO; Supplementary Table S1) with 1,878 patients of European ancestry. To increase statistical power, Set 2 was combined with Set 1 via meta-analysis (due to separate genotyping). Patients with incomplete survival time information (e.g., missing vital status, enrollment greater than 10 years after diagnosis) were excluded, resulting in 1,744 patients available for analysis (1,027 deaths at 10 years postdiagnosis). Sixteen patients had missing follow-up time, and were set to the median follow-up time of the study; 6 patients had incomplete time-to-entry information, and were set to entry at time zero. Characteristics of both sample sets are described in Table 1. Both patient sets consist of predominantly patients with high-grade serous histology, although Set 1 also includes a large number of endometrioid patients.
. | Sample Set 1 . | Sample Set 2 . | ||||
---|---|---|---|---|---|---|
. | Alive (N = 2,036) . | Deceased (N = 2,257) . | Total (N = 4,293) . | Alive (N = 717) . | Deceased (N = 1,027) . | Total (N = 1,744) . |
Median survival time (years, 95% CI) | 4.91 (4.65–5.21) | 2.56 (2.44–2.68) | 4.58 (4.41–4.77) | 6.03 (5.70–6.43) | 2.89 (2.69–3.08) | 4.32 (4.10–4.60) |
Person years (years) | 10,737.46 | 6,169.05 | 16,906.51 | 3,545.56 | 2,805.42 | 6,350.98 |
Grade | ||||||
Well differentiated | 288 (14.1%) | 79 (3.5%) | 367 (8.5%) | 69 (9.6%) | 72 (7.0%) | 141 (8.1%) |
Moderately differentiated | 426 (20.9%) | 370 (16.4%) | 796 (18.5%) | 109 (15.2%) | 216 (21.0%) | 325 (18.6%) |
Poorly differentiated | 886 (43.5%) | 1,318 (58.4%) | 2,204 (51.3%) | 3,237 (47.0%) | 554 (53.9%) | 891 (51.1%) |
Undifferentiateda | 276 (13.6%) | 302 (13.4%) | 578 (13.5%) | 0 (0.0%) | 0 (0.0%) | 0 (0.0%) |
Unknown | 160 (7.9%) | 188 (8.3%) | 348 (8.1%) | 202 (28.2%) | 185 (18.0%) | 387 (22.2%) |
Stage | ||||||
Localized | 436 (21.4%) | 88 (3.9%) | 524 (12.2%) | 111 (15.5%) | 32 (4.5%) | 143 (8.2%) |
Regional | 444 (21.8%) | 163 (7.2%) | 607 (14.1%) | 216 (30.1%) | 239 (33.3%) | 155 (26.1%) |
Distant | 1,091 (53.6%) | 1,949 (86.4%) | 3,040 (70.8%) | 244 (34.0%) | 635 (88.6%) | 879 (50.4%) |
Unstaged/unknown | 65 (3.2%) | 57 (2.5%) | 122 (2.8%) | 146 (20.4%) | 121 (16.9%) | 267 (15.3%) |
Histology | ||||||
Serous | 1,349 (66.3%) | 1,943 (86.1%) | 3,292 (76.7%) | 706 (98.5%) | 1,018 (99.1%) | 1,724 (98.9%) |
Endometrioid | 536 (26.3%) | 199 (8.8%) | 735 (17.1%) | 2 (0.3%) | 2 (0.2%) | 4 (0.2%) |
Mucinous | 37 (1.8%) | 5 (0.2%) | 42 (1.0%) | 0 (0.0%) | 1 (0.1%) | 1 (0.1%) |
Clear cell | 25 (1.2%) | 8 (0.4%) | 33 (0.8%) | 8 (1.1%) | 3 (0.3%) | 11 (0.6%) |
Mixed cell | 16 (0.8%) | 21 (0.9%) | 37 (0.9%) | 1 (0.1%) | 2 (0.2%) | 3 (0.2%) |
Other epithelial | 73 (3.6%) | 81 (3.6%) | 154 (3.6%) | 0 (0.0%) | 1 (0.1%) | 1 (0.1%) |
Age, mean (SD) | 56.8 (11.0) | 60.7 (10.8) | 58.8 (11.1) | 56.7 (10.6) | 60.2 (9.7) | 58.8 (10.2) |
. | Sample Set 1 . | Sample Set 2 . | ||||
---|---|---|---|---|---|---|
. | Alive (N = 2,036) . | Deceased (N = 2,257) . | Total (N = 4,293) . | Alive (N = 717) . | Deceased (N = 1,027) . | Total (N = 1,744) . |
Median survival time (years, 95% CI) | 4.91 (4.65–5.21) | 2.56 (2.44–2.68) | 4.58 (4.41–4.77) | 6.03 (5.70–6.43) | 2.89 (2.69–3.08) | 4.32 (4.10–4.60) |
Person years (years) | 10,737.46 | 6,169.05 | 16,906.51 | 3,545.56 | 2,805.42 | 6,350.98 |
Grade | ||||||
Well differentiated | 288 (14.1%) | 79 (3.5%) | 367 (8.5%) | 69 (9.6%) | 72 (7.0%) | 141 (8.1%) |
Moderately differentiated | 426 (20.9%) | 370 (16.4%) | 796 (18.5%) | 109 (15.2%) | 216 (21.0%) | 325 (18.6%) |
Poorly differentiated | 886 (43.5%) | 1,318 (58.4%) | 2,204 (51.3%) | 3,237 (47.0%) | 554 (53.9%) | 891 (51.1%) |
Undifferentiateda | 276 (13.6%) | 302 (13.4%) | 578 (13.5%) | 0 (0.0%) | 0 (0.0%) | 0 (0.0%) |
Unknown | 160 (7.9%) | 188 (8.3%) | 348 (8.1%) | 202 (28.2%) | 185 (18.0%) | 387 (22.2%) |
Stage | ||||||
Localized | 436 (21.4%) | 88 (3.9%) | 524 (12.2%) | 111 (15.5%) | 32 (4.5%) | 143 (8.2%) |
Regional | 444 (21.8%) | 163 (7.2%) | 607 (14.1%) | 216 (30.1%) | 239 (33.3%) | 155 (26.1%) |
Distant | 1,091 (53.6%) | 1,949 (86.4%) | 3,040 (70.8%) | 244 (34.0%) | 635 (88.6%) | 879 (50.4%) |
Unstaged/unknown | 65 (3.2%) | 57 (2.5%) | 122 (2.8%) | 146 (20.4%) | 121 (16.9%) | 267 (15.3%) |
Histology | ||||||
Serous | 1,349 (66.3%) | 1,943 (86.1%) | 3,292 (76.7%) | 706 (98.5%) | 1,018 (99.1%) | 1,724 (98.9%) |
Endometrioid | 536 (26.3%) | 199 (8.8%) | 735 (17.1%) | 2 (0.3%) | 2 (0.2%) | 4 (0.2%) |
Mucinous | 37 (1.8%) | 5 (0.2%) | 42 (1.0%) | 0 (0.0%) | 1 (0.1%) | 1 (0.1%) |
Clear cell | 25 (1.2%) | 8 (0.4%) | 33 (0.8%) | 8 (1.1%) | 3 (0.3%) | 11 (0.6%) |
Mixed cell | 16 (0.8%) | 21 (0.9%) | 37 (0.9%) | 1 (0.1%) | 2 (0.2%) | 3 (0.2%) |
Other epithelial | 73 (3.6%) | 81 (3.6%) | 154 (3.6%) | 0 (0.0%) | 1 (0.1%) | 1 (0.1%) |
Age, mean (SD) | 56.8 (11.0) | 60.7 (10.8) | 58.8 (11.1) | 56.7 (10.6) | 60.2 (9.7) | 58.8 (10.2) |
aGrade was based on a 4-tier grading system for Set 1 and a 3-tier grading system for Set 2.
Genotyping and quality control
For Set 1, genotyping was performed using genomic DNA and whole genome–amplified (WGA) DNA for 5,904 and 389 patients, respectively. Genotyping was performed at Affymetrix Corporation using the Affymetrix Axiom Exome Array (www.affymetrix.com) of 409,582 markers, including standard content markers, as well as 100,000 custom markers chosen based on preliminary survival associations, follow-up of prior risk associations, and noncoding candidate gene variants. The default quality control (QC) criteria recommended by Affymetrix was applied to all samples, and then the WGA and genomic samples were processed separately. For genomic samples, genotypes were recalled using Powertool for variants with MAF < 5%, and samples and variants with call rates < 97% were excluded. WGA samples and variants with call rates < 97% were also excluded and then merged with genomic samples, resulting in 384,029 variants. Excluding monomorphic variants, 227,892 variants were used in analysis.
For Set 2, patients contributed genomic DNA samples that were genotyped at the University of Cambridge (Cambridge, United Kingdom) on the Illumina Infinium HumanExome BeadChip (www.Illumina.com) at a total of 247,840 exonic markers. Genotype calling was carried out according to Best Practice Guidelines using the GenCall module in Illumina's Genome Studio with a default GenCall threshold of 0.15 (25). Samples with low call rate (<99%), high or low heterozygosity based on common SNPs (MAF ≥ 0.05), and ambiguous sex or relatedness were excluded. A total of 130,909 variants were monomorphic and 147 variants were excluded due to low call rate, resulting in 114,620 markers.
Of the 73,203 variants overlapping on both genotyping arrays, a majority are rare variants (median MAF = 0.0028; IQR = 0.068) and coding variants (89.2%). Supplementary Fig. S1 displays the MAF distribution by functional category.
Single variant analysis
Both patient sets were analyzed using Cox proportional hazards regression to estimate the association of each variant with overall 10-year survival time. For each variant, a Cox proportional hazards model was fit to the number of copies of the minor allele, and a likelihood ratio test was conducted, accounting for left truncation (26) and right-censoring after 10 years. For top variants, Kaplan–Meier curves were also generated that account for left truncation and right censoring. Because the primary and secondary patient sets were genotyped separately and on different platforms, they were analyzed separately.
For the primary patient set, models were first minimally adjusted for age, site, and three principal components (PC) to account for population stratification. PCs were computed on the basis of common variants (MAF > 1%, HWE P > 1E−7) using the study patients, as well as HapMap subjects. In addition, the following clinical factors with univariate survival associations of P < 0.05 and with P < 0.0001 in multivariable modeling were included as covariates along with three PCs: age, stage, grade, histology, and site. Both sample sets consisted of patients of primarily serous histology, although Set 1 also included many endometrioid cases. Therefore all histologies were included in analyses (N = 4293), but histologic subtype-specific analyses were also conducted for high-grade serous (N = 3149) and endometrioid subtypes (N = 735), because some candidate variants were targeted for these subtypes. Variants with P < 2.2E−7 were considered significant after Bonferroni correction for multiple testing.
Because the sample size of Set 2 was relatively small, its primary utility was to improve power and provide supportive evidence for the primary results via meta-analysis; the patient sets were genotyped separately, which limited the ability to conduct a pooled analysis. For Set 2, Cox proportional hazards analyses (accounting for left truncation and right-censoring after 10 years) were conducted using minimal adjustments as above; no other covariates were associated with the outcome and therefore additional adjustments were not included due to the smaller sample size. Fixed effect meta-analyses were conducted for the overlapping variants on both the Affymetrix and Illumina platforms (73,203 variants), based on the likelihood ratio summary statistics from the minimally adjusted Set 1 and Set 2 analyses. Single variant meta-analysis was conducted using the R package “rMeta” (http://cran.r-project.org/web/packages/rmeta/). All analyses were conducted in R version 3.0.2 (http://www.R-project.org/).
Gene-level analysis
As single variant analyses are under-powered for rare variants, gene-level tests were also used. Variants were mapped to genes based on Human Genome Build 37 and annotation supplied by Affymetrix and Illumina, resulting in 18,323 genes analyzed in Set 1 and 13,191 genes analyzed in Set 2. To be included in gene-level analysis, genes were required to contain at least two SNPs. Gene-level testing was conducted on both the sample sets separately using burden tests, as well as Sequence Kernel Association Test (SKAT) for Cox proportional hazards models (27). The gene-level burden test is a weighted sum of the genotypes of all variants with a gene regressed on survival time (20). The gene-level SKAT statistic is a weighted sum of the single variant likelihood ratio test statistics, across all variants within a gene (23). Variant weights were computed as a function of MAF, where all variants were included, and rare variants received much greater weight; specifically, the weight function used was |$\sqrt {w_j} = Beta\left({MAF_j, \, 1,\,25} \right)$| for each variant j (23). Minimally adjusted (Sets 1 and 2) and fully adjusted tests (Set 1) were computed as described above for single variant tests. Gene-level tests (SKAT and burden tests) were computed using the R package “seqMeta” (http://cran.r-project.org/web/packages/seqMeta/), with minor modifications to accommodate left-truncated data.
As with the single variant tests, meta-analyses of gene-level tests across the sample sets were also conducted. For the purpose of meta-analysis only, variants that were monomorphic in one set but not the other were included to provide information on MAF across both studies. A total of 13,163 genes available in both sample sets were meta-analyzed using the R-package “seqMeta”. Genes with P < 3.8E−6 were considered significant after Bonferroni correction for multiple testing. Sensitivity analyses were conducted for variants with MAF < 0.05, MAF < 0.10, and MAF < 0.25, as well as for coding variants only.
Results
Single variants
Results for all 73,203 variants included in the meta-analysis of 6,054 invasive EOC patients are displayed in Fig. 1. Three variants were associated with overall survival at P values <1.0E−4, including two with MAF < 0.001 (Table 2). The top ranking variant in the meta-analysis was common variant rs8170 in BABAM1 at 19p13.11 (Pmeta = 1.14E−6), previously reported to be associated with EOC risk (10); the minor allele of this variant was associated with decreased survival in both sets (HR1 = 1.17, HR2 = 1.14; Fig. 2). Additional top ranking variants were rare nonsynonymous variants rs140079492 in PAK7 (P21 protein (cdc42/rac)-activated kinase) at 20p12.2 [(Pmeta = 9.32E−6) and rs34335714 in CTSW (cathepsin W) at 11q13.1 (Pmeta = 4.91E−5, Table 2)].
. | . | . | . | Set 1 (N = 4,293) . | . | . | Set 2 (N = 1,744) . | . | . | Meta-analysis . | |||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
Chr . | rsID (allele) . | Gene . | Annotation . | HR (95% CI) . | P . | MAF . | Total number of patients with minor allele . | Number of deceased patients with minor allele . | HR (95% CI) . | P . | MAF . | Total number of patientswith minor allele . | Number of deceased patients with minor allele . | HR (95% CI) . | Meta P . |
11 | rs34335714 (A/G) | CTSW | Nonsynonymous | 5.34 (2.00–14.25) | 0.0008 | 0.0007 | 6 | 4 | 10.81 (1.52–76.79) | 0.0174 | 0.0003 | 1 | 1 | 6.15 (2.56–14.79) | 4.91E–5 |
19 | rs8170 (A/G) | BABAM1 | Previous GWAS hit | 1.17 (1.09–1.26) | 2.23E–5 | 0.1966 | 1446 | 839 | 1.14 (1.03–1.26) | 0.0151 | 0.2124 | 664 | 433 | 1.16 (1.09–1.23) | 1.14E–6 |
20 | rs140079492 (C/T) | PAK7 | Nonsynonymous | 3.41 (1.42–8.19) | 0.0062 | 0.0006 | 5 | 5 | 8.52 (2.74–26.45) | 0.0002 | 0.0011 | 4 | 3 | 4.80 (2.40–9.61) | 9.32E–6 |
. | . | . | . | Set 1 (N = 4,293) . | . | . | Set 2 (N = 1,744) . | . | . | Meta-analysis . | |||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
Chr . | rsID (allele) . | Gene . | Annotation . | HR (95% CI) . | P . | MAF . | Total number of patients with minor allele . | Number of deceased patients with minor allele . | HR (95% CI) . | P . | MAF . | Total number of patientswith minor allele . | Number of deceased patients with minor allele . | HR (95% CI) . | Meta P . |
11 | rs34335714 (A/G) | CTSW | Nonsynonymous | 5.34 (2.00–14.25) | 0.0008 | 0.0007 | 6 | 4 | 10.81 (1.52–76.79) | 0.0174 | 0.0003 | 1 | 1 | 6.15 (2.56–14.79) | 4.91E–5 |
19 | rs8170 (A/G) | BABAM1 | Previous GWAS hit | 1.17 (1.09–1.26) | 2.23E–5 | 0.1966 | 1446 | 839 | 1.14 (1.03–1.26) | 0.0151 | 0.2124 | 664 | 433 | 1.16 (1.09–1.23) | 1.14E–6 |
20 | rs140079492 (C/T) | PAK7 | Nonsynonymous | 3.41 (1.42–8.19) | 0.0062 | 0.0006 | 5 | 5 | 8.52 (2.74–26.45) | 0.0002 | 0.0011 | 4 | 3 | 4.80 (2.40–9.61) | 9.32E–6 |
NOTE: Analyses are adjusted for age, site, and three PCs.
Despite reduced sample size compared with meta-analysis, Set 1–specific results are also of interest due to the much larger number of variants targeted than in Set 2 (including custom candidate variants); therefore, we also report the results of all variants analyzed in Set 1 (including those non-overlapping with Set 2). In this primary EOC sample set (4,293 patients), with adjustment for age, site, stage, grade, histology, and PCs, 24 variants in 11 chromosomal regions showed evidence of association with overall survival with P < 5.0E−5, although none of these associations are significant after correction for multiple testing. Genome-wide results are plotted in Supplementary Fig. S2 and top variants are described in Supplementary Table S2. The variant most strongly associated with survival was rs7642051, a common variant (MAF = 0.46) in LMCD1-AS1 (LMCD1 antisense RNA 1) at 3p26.1 (P = 4.10E−6); its minor allele was associated with decreased survival (HR = 1.15; Supplementary Table S2, Supplementary Fig. S3). Several correlated variants in an intergenic region of chromosome 1 between UBE2U and CACHD1 (1p31.3), targeted due to a prior study of overall survival, also showed modest evidence for association with survival (Supplementary Fig. S4). This region showed similar results for the high-grade serous subtype (data not shown). Among patients with endometrioid subtype, the most strongly associated variant was rs757759 in TBXAS1 at 7q34 (HR = 1.78, P = 1.55E−6, MAF = 0.23).
Gene-level
Results for all 13,163 genes included in the meta-analysis across Sets 1 and 2 are displayed in Fig. 3. Five genes had meta-analysis P < 1.0E−4 based on the burden test, and four of these genes showed evidence of association in both sets separately (Table 3). The top-ranking gene was ATG2B (Pmeta = 1.08E−6), which was significant after multiple testing correction based on 13,163 genes (Pcorrected = 0.014). Out of 17 single variants meta-analyzed in ATG2B at 14q32.2, nine had P < 0.10 (Supplementary Fig. S5) and 15 with positive effect size estimates (HR > 1.0); of note, variants were largely uncorrelated (pairwise r2 < 0.2), and most were rare (16 with MAF < 0.01; median = 5.9E−4, range = 1.7E−4–6.8E−3). When restricted to the high-grade serous subset, ATG2B remains the top gene (Pmeta = 1.75E−6). Other genes of interest across all histologies include zinc finger protein PEG3 (Pmeta = 1.82E−5), helicase-like transcription factor HLTF (Pmeta = 1.9E−5), mitochondrial fission regulator 1-like gene MTFR1L (P = 3.1E−5), and T-cell surface glycoprotein CD1E (Pmeta = 8.78E−5; Table 3). Results differed on the basis of the SKAT meta-analysis; while gene rankings were relatively similar to the burden test, P values were higher (min P = 1.17E−4 for ZNF131). In particular, the meta-analysis P value based on SKAT for ATGB2 was P = 0.006 (Table 3). Although top-ranked in single variant analyses, BABAM1 was not significant at the gene-level (Pburden = 0.83, PSKAT = 0.53).
. | . | Set 1 (N = 4,293) . | . | Set 2 (N = 1,744) . | . | Burden . | SKAT . | ||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|
Chr . | Gene . | N variants . | CMAF . | Median MAF . | N variants . | CMAF . | Median MAF . | P Set 1 . | P Set 2 . | Pmeta . | P Set 1 . | P Set 2 . | Pmeta . |
1 | MTFR1L | 8 | 0.0032 | 0.00024 | 5 | 0.0026 | 0.00057 | 2.21E–03 | 0.0028 | 3.13E–05 | 0.1089 | 0.0035 | 0.0037 |
1 | CD1E | 12 | 0.5692 | 0.00047 | 3 | 0.3374 | 0.00057 | 9.29E–05 | 0.7255 | 8.78E–05 | 0.0002 | 0.7184 | 0.0002 |
3 | HLTF | 17 | 0.7604 | 0.00128 | 7 | 0.0801 | 0.00659 | 1.09E–03 | 0.0049 | 1.91E–05 | 0.0310 | 0.0538 | 0.0096 |
14 | ATG2B | 34 | 1.7512 | 0.00058 | 17 | 0.4765 | 0.00057 | 8.63E–05 | 0.0030 | 1.08E–06 | 0.0790 | 0.0078 | 0.0057 |
19 | PEG3 | 38 | 1.0465 | 0.00058 | 18 | 0.2036 | 0.00143 | 3.18E–04 | 0.0221 | 1.82E–05 | 0.0630 | 0.0856 | 0.0751 |
. | . | Set 1 (N = 4,293) . | . | Set 2 (N = 1,744) . | . | Burden . | SKAT . | ||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|
Chr . | Gene . | N variants . | CMAF . | Median MAF . | N variants . | CMAF . | Median MAF . | P Set 1 . | P Set 2 . | Pmeta . | P Set 1 . | P Set 2 . | Pmeta . |
1 | MTFR1L | 8 | 0.0032 | 0.00024 | 5 | 0.0026 | 0.00057 | 2.21E–03 | 0.0028 | 3.13E–05 | 0.1089 | 0.0035 | 0.0037 |
1 | CD1E | 12 | 0.5692 | 0.00047 | 3 | 0.3374 | 0.00057 | 9.29E–05 | 0.7255 | 8.78E–05 | 0.0002 | 0.7184 | 0.0002 |
3 | HLTF | 17 | 0.7604 | 0.00128 | 7 | 0.0801 | 0.00659 | 1.09E–03 | 0.0049 | 1.91E–05 | 0.0310 | 0.0538 | 0.0096 |
14 | ATG2B | 34 | 1.7512 | 0.00058 | 17 | 0.4765 | 0.00057 | 8.63E–05 | 0.0030 | 1.08E–06 | 0.0790 | 0.0078 | 0.0057 |
19 | PEG3 | 38 | 1.0465 | 0.00058 | 18 | 0.2036 | 0.00143 | 3.18E–04 | 0.0221 | 1.82E–05 | 0.0630 | 0.0856 | 0.0751 |
NOTE: Analyses are adjusted for age, site, and three PCs.
Abbreviation: CMAF, combined (sum) minor allele frequency across all variants in the gene.
In analysis of Set 1 alone, when adjusted for age, site, stage, grade, histology, and PCs, five genes showed evidence of association with overall survival at P < 1.0 × 10−4 from the burden test, although none of these associations are significant after multiple testing correction based on number of genes (Supplementary Table S3); results were similar using SKAT (Supplementary Fig. S6). The top-ranked genes were POGLUT1 (protein O-glucosyltransferase 1; Pburden = 7.10E−6, PSKAT = 8.33E−6), ST20 (suppressor of tumorigenicity 20; Pburden = 2.32E−5, PSKAT = 2.32E−5), and ATG2B (autophagy-related 2B; Pburden = 5.16E−5, PSKAT = 3.23E−2), the most significant gene in meta-analysis. In addition, consistent with individual variant analysis, variation in UBE2U showed some evidence of association (36 variants, Pburden = 5.07E−4, PSKAT = 4.95E−4), although evidence was less strong in the CACHD1 region (65 variants, Pburden = 0.004, PSKAT = 0.025). Gene-level results were similar when restricted to patients with high-grade serous subtype, when restricted to variants with MAF < 0.05, MAF < 0.10, and MAF < 0.25, and when restricted to only coding variants (data not shown); ATG2B remained the top gene, with meta-analysis P values based on a burden test ranging from 3.94E−7 to 1.08E−6.
Discussion
To our knowledge, this study represents the first exome-wide assessment of association between EOC survival and rare genetic variation, both for single variants and for genes, including over 6,000 patients. By combining data from two exome-based OCAC genotyping projects, genome-wide significant findings were identified at the gene-level in ATG2B (Pmeta = 1.1E−6, Pcorrected = 0.01); and, in analysis of approximately 73,000 variants, the known EOC susceptibility variant rs8170 in BABAM1 arose as the most statistically significant (Pmeta = 1.1E−6), followed by nonsynonymous variants rs140079492 in PAK7 (serine/threonine-protein kinase; missense mutation Glu → Gly) and rs34335714 in CTSW (cathepsin W; missense mutation Ser → Asn, predicted splice variant). As genotype imputation methods are known to be unreliable in the context of rare variation (28, 29), the use of direct genotyping is an important strength of this report.
ATG2B, the survival-associated gene based on burden testing, is a key autophagy gene residing on chromosome 14q32.2. Autophagy is a cell death process which uses degradation of lysosomal cytoplasmic components; in cancer, it is thought to link apoptosis with programmed necrosis and has been proposed as an alternative target to treat tumor resistance (30). It is regulated by several miRNAs, including mir-30d and mir-130a (31, 32), and has recently been linked to precursors of pancreatic cancer (33). Here, up to 34 ATG2B variants genotyped in either Sets 1 or 2 associated in combination with EOC survival; this was consistent with the trends of individual ATG2B variants, where 9 of 17 ATG2B variants observed in both sets had P < 0.10 in the meta-analysis. These results suggest that ATG2B may warrant targeted sequencing in large datasets to confirm the existence of potentially rare, inherited prognostic variants.
Although not statistically significant after correction for multiple testing, BABAM1 rs8170 at 19p13, the top ranking individual variant across both sample sets (adjusting for age, site, and PCs), warrants attention. This synonymous coding variant is relatively common (MAF = 0.20) and is included on multiple commercial GWAS arrays, as well as the Affymetrix and Illumina exome-based arrays used here. Its MAF contributed to increased power to its detection over the rarer variants on the arrays. BABAM1 rs8170 is known to associate with EOC risk (10), with triple-negative breast cancer risk (34), and with risk of ovarian cancer and estrogen-negative subtypes of breast cancer among BRCA1 and BRCA2 mutation carriers (16, 35). The estimated hazard ratio in meta-analysis was 1.16, which roughly translates to a 16% reduction in median survival time for patients with one copy of the minor allele compared with those with no copies (36). Sensitivity analysis of the current data show that when also adjusted for stage, grade, and histology, the rs8170 effect size was slightly reduced (HR = 1.11 vs. HR = 1.17; P = 4.33E−5). However, when the analysis was restricted to patients with high-grade serous cancer, effect size estimates were relatively robust whether adjusted for age, site, and PCs (HR = 1.17, P = 4.66E−5) or also stage and grade (HR = 1.14, P = 8.02E−5). Additional patient subset analyses may uncover other patterns of survival; for example, for a subset of OCAC with available data regarding detailed treatment strategies, no association of rs8170 with overall survival or progression-free survival was observed after stratifying by treatment (37). Although prior EOC survival studies yielded inconsistent results (10), unpublished meta-analysis of genotyped and imputed variation at rs8170 from a much larger sample of over 18,000 EOC patients (including the majority of those in the current study and all of the patients from a prior EOC survival study; ref. 10) revealed a genome-wide significant association (HR = 1.15, P = 4.7E−9; Pharoah and colleagues, personal communication). Altogether, the BABAM1 region containing rs8170 is clearly of importance to ovarian cancer and will be subjected to fine-mapping and detailed functional analyses.
As one of the early GWAS using exome-based Affymetrix and Illumina genotyping arrays, we note some rare variant analytical lessons learned. Single variant analyses often yielded different results than the aggregate gene-level tests, and represent useful, complementary analysis strategies. Because single variant tests were subject to reduced power to detect associations with rare variants, only associations with common variants (such as rs8170) and rare variants with extremely large effect sizes were detectable. For example, in the single variant meta-analysis, the median MAF is less than 0.005, which translates to less than 60 subjects out of the combined Set 1 and Set 2 sample size of 6,037 that can be expected to harbor at least one copy of the rare allele. For rare variants with MAF of 0.005, we have 80% power to detect HRs greater than 2.9, and for even rarer variants, only much larger HRs would be detectable with this sample size. Thus, larger sample sets will be necessary to rule out or detect individual associations with the approximately 70,000 variants that were individually meta-analyzed here. However, because gene-level tests consider all variants within a region and weight rare variants more heavily, they are able to identify associations with regions that contain many variants with moderate evidence of association; for example, no variant in ATG2B had a meta-analysis P < 0.001, yet it was genome-wide significant at the gene-level.
The most powerful gene-level analysis tool is dependent on the underlying genetic etiology of the trait, and because the true genetic architecture of EOC survival is unknown, we used two methods for gene-level testing that have been shown to be powerful under different scenarios (20, 27). While the gene rankings using a burden test or SKAT were similar, the P values based on the SKAT analysis were often higher, suggesting that SKAT may be less powerful, consistent with simulation studies (23). The burden test assumes that the direction of effect for each variant within the gene region is the same, while SKAT does not assume consistent effect direction. Therefore the burden test will be more powerful to detect rare variants if the true causal effects are in the same direction (23); because our outcome of interest in this study was survival, it is not unreasonable to assume that rare variants would increase risk of death (rather than the converse), although of course the true underlying genetic model of EOC survival remains unknown. In fact, for ATG2B (where the SKAT P value was higher than the burden P value), all but one of the nine variants with Pmeta < 0.10 had estimated HRs greater than 1.0.
The findings of this study should be considered in the context of the following limitations, in addition to sample size. The two patient sets were genotyped separately on different platforms that targeted different variants due to differences in chemistry (Affymetrix vs. Illumina), even though both panels were designed from a common set of variants from the Exome Sequencing Project (24). This is far from ideal, not only requiring extensive additional data harmonization measures, but limiting our ability to conduct a pooled analysis and restricting meta-analysis to a reduced set of variants (only 73,203 combined from 227,892 Set 1 and 114,620 Set 2 variants) and genes (only 13,163 combined from 18,323 Set 1 and 13,191 from Set 2 genes) that were observed in both sets. Therefore many rare variants that were observed in the larger Set 1 were not assessed in Set 2. Furthermore, although also a strength, approximately 25% of the Affymetrix array was customized content, and therefore not interrogated in Set 2, including the intergenic region in chromosome 1 between UBE2U and CACHD1. Imputation to a denser set of variants, for example using 1000 Genomes Project data, is known to be unreliable based on rare variants and was therefore not performed; future work to integrate these data with other available genotypes for these patients will be informative. In addition, while all patients underwent surgery, the important clinical factors of optimal debulking and type of treatment (e.g., chemotherapy) could not be included as covariates or stratification factors due to a large amount of missing data; both should be examined in future genetic studies to assess their impact on survival outcomes. Finally, it should be noted that some of the data analyzed here were not independent from previous OCAC EOC overall survival reports (5, 10, 38); specifically, some variants on the exome arrays were also on GWAS or candidate gene arrays which had previously been used on a subset of the current patients. However, the vast majority of the current array content was largely used to test the novel rare variant hypothesis.
Finally, this study comprised data from 20 OCAC studies, with differing ascertainment strategies and methods for clinical follow-up of cases, resulting in the potential for increased study heterogeneity in the overall survival outcome measure. While we adjusted for study site in all models, uncaptured heterogeneity across studies (such as differences in treatment strategies over time and across medical institutions) will lead to increased noise and reduced power. In future studies, more refined outcome measures targeted at reducing this heterogeneity could improve power. For example, because first-line therapy may differ from treatment strategies after a recurrence, progression-free survival may be a more optimal endpoint, although this measure is also susceptible to study site heterogeneity. Because recurrence data was not widely available for all patients across all sites, it was not considered as a primary endpoint in this study. However, analysis of the association between rs8170 and progression-free survival were consistent with the results observed for the overall survival outcome (HRmeta = 1.15, Pmeta = 1.3E−6, 95% CI = 1.09–1.22).
To summarize, in the absence of associations between common genetic variants and EOC survival, we examined the role of primarily rare exonic variants, individually and collectively across genes, using partially customized commercial arrays in over 6,000 patients. While no individual variants were significant at the genome-wide level, ATG2B showed gene-level evidence of association with overall survival, and BABAM1 rs8170 was again highlighted to be of particular relevance. This suggests potential candidates for future studies that may lead to targets for improved EOC outcomes, although follow-up in the larger sample of OCAC subjects not included here, as well as a focus on reducing the heterogeneity in the outcome measure, will be critical.
Disclosure of Potential Conflicts of Interest
No potential conflicts of interest were disclosed.
Authors' Contributions
Conception and design: S.J. Winham, Y.A. Chen, H. Anton-Culver, H.-Y. Lin, R.A. Vierkant, A. Ziogas, A. Berchuck, P.D.P. Pharoah, B.L. Fridley, T.A. Sellers
Development of methodology: S.J. Winham, H. Anton-Culver
Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): M.C. Larson, H. Anton-Culver, E.V. Bandera, D. Cramer, J.A. Doherty, M.T. Goodman, J. Gronwald, B.Y. Karlan, S.K. Kjaer, D.A. Levine, U. Menon, R.B. Ness, T. Pejovic, M. A. Rossing, N. Wentzensen, Y.T. Bean, M.E. Carney, J.M. Cunningham, C. Cybulski, A. deFazio, R.P. Edwards, S.A. Gayther, A. Gentry-Maharaj, M. Gore, A. Jensen, J. Lester, J. Lissowska, J. Lubinski, J. Menkiszak, F. Modugno, K.B. Moysich, I. Orlow, M.C. Pike, S.J. Ramus, H. Song, K.L. Terry, P.J. Thompson, D.J. van den Berg, A.F. Vitonis, C. Walsh, A.H. Wu, A. Ziogas, G. Chenevix-Trench, J.M. Schildkraut, C.M. Phelan, P.D.P. Pharoah, T.A. Sellers
Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): S.J. Winham, A. Pirie, Y.A. Chenm, Z.C. Fogarty, N. Wentzensen, E.M. Dicks, E.S. Iversen, J. Menkiszak, K.B. Moysich, J.P. Tyrer, H. Yang, A. Ziogas, P.D.P. Pharoah, B.L. Fridley
Writing, review, and/or revision of the manuscript: S.J. Winham, A. Pirie, Y.A. Chen, M.A. Earp, H. Anton-Culver, E.V. Bandera, D. Cramer, J.A. Dohertym, M.T. Goodman, J. Gronwald, B.Y. Karlan, S.K. Kjaer, D.A. Levine, U. Menon, R.B. Ness, C.L. Pearce, M.A. Rossing, N. Wentzensen, L.A. Brinton, M.E. Carney, J.M. Cunningham, C. Cybulski, A. deFazio, R.P. Edwards, S.A. Gayther, A. Gentry-Maharaj, M. Gore, A. Jensen, H.-Y. Lin, J. Lissowska, F. Modugno, K.B. Moysich, I. Orlow, S.J. Ramus, H. Song, R.A. Vierkant, C. Walsh, A.H. Wu, H. Yang, A. Berchuck, G. Chenevix-Trench, J. Permuth-Wey, C.M. Phelan, P.D.P. Pharoah, B.L. Fridley, T.A. Sellers, E.L. Goode
Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): M.C. Larson, B.Y. Karlan, S.E. Johnatty, J. Lissowska, J. Lubinski, F. Modugno, S.J. Ramus, P.J. Thompson, A.F. Vitonis, H. Yang, A. Ziogas, A. Berchuck
Study supervision: S.K. Kjaer, E.S. Iversen, J. Lissowska, F. Modugno, K.B. Moysich, P.J. Thompson, B.L. Fridley
Other (preparation and submitting of DNA samples): M. Bisogna
Other (principal investigator of primary source of funding grant): T.A. Sellers
Acknowledgments
The Australian Ovarian Cancer Study Management Group (D. Bowtell, G. Chenevix-Trench, A. deFazio, D. Gertig, A. Green, P. Webb) and ACS Investigators (A. Green, P. Parsons, N. Hayward, P. Webb, D. Whiteman) thank all the clinical and scientific collaborators (see http://www.aocstudy.org/) and the women for their contribution. The authors thank I. Jacobs, M. Widschwendter, E. Wozniak, A. Ryan, J. Ford, and N. Balogun for their contribution to the study.
Grant Support
This work was supported by grants from the NIH Office of Research on Women's Health (Building Interdisciplinary Careers in Women's Health award K12HD065987; to S.J. Winham) and Genetic Associations and Mechanisms in Oncology (GAME-ON), an NCI Cancer Post-GWAS Initiative (U19-CA148112; to T.A. Sellers). G. Chenevix-Trench is supported by Fellowships from NHMRC. A. deFazio is supported by the University of Sydney Cancer Research Fund and the Cancer Institute NSW through the Sydney West Translational Cancer Research Centre. C.M. Phelan is funded in part by NIH R01 CA-149429. A. Pirie is funding in part by a Medical Research Council Studentship. A. Gentry-Maharaj is supported by Eve Appeal and the Oak Foundation. In addition, the individual study sites were supported by: AUS: U.S. Army Medical Research and Materiel Command (DAMD17-01-1-0729; to G. Chenevix-Trench), National Health & Medical Research Council of Australia (199600 and 400281; to G. Chenevix-Trench), Cancer Councils of New South Wales, Victoria, Queensland, South Australia, and Tasmania, Cancer Foundation of Western Australia. DOV: Funding for this study was provided by U19-CA148112 (to T. Sellers); HAW: NIH (R01-CA58598, N01-CN-55424, and N01-PC-67001; to M.T. Goodman); HOP: US Army Medical Research and Material Command DAMD17-02-1-0669 (to F. Modugno); NCI K07-CA080668, R01-CA95023, R01-CA126841 (to F. Modugno); P50-CA159981 (to F. Modugno); NIH/National Center for Research Resources/General Clinical Research Center grant M01-RR000056 (to F. Modugno); LAX: American Cancer Society Early Detection Professorship (SIOP-06-258-01-COUN; to B.Y. Karlan) and the National Center for Advancing Translational Sciences (NCATS), grant UL1TR000124 (to B.Y. Karlan); MAL: funding for this study was provided by research grant R01- CA61107 from the NCI, Bethesda, MD (to S.K. Kjaer and A. Jensen); research grant 94 222 52 from the Danish Cancer Society, Copenhagen, Denmark (to S.K. Kjaer and A. Jensen); and the Mermaid I project (to S.K. Kjaer and A. Jenesn); MAC and MAY: NIH (R01-CA122443, P30-CA15083, P50-CA136393; to E.L. Goode); NCO: Department of Defense (DAMD17-02-1-0666; to J.M. Schildkraut); NEC: NIH R01-CA54419 and P50-CA105009 and Department of Defense W81XWH-10-1-02802 (to K.L. Terry); NJO: NCI (NIH-K07 CA095666, NIH-K22-CA138563, and P30-CA072720; to E.V. Bandera) and the Cancer Institute of New Jersey (to E.V. Bandera); ORE: OHSU Foundation (to T. Pejovic); POC: Pomeranian Medical University (to J. Gronwald); RMH: Cancer Research UK (no grant number is available), Royal Marsden Hospital (to P.D. Pharoah); SEA: Cancer Research UK (C490/A10119 C490/A10124; to P.D. Pharoah); UK National Institute for Health Research Biomedical Research Centres at the University of Cambridge, SEARCH team, Craig Luccarini, Caroline Baynes, Don Conroy; UKO: The UKOPS study was funded by The Eve Appeal (The Oak Foundation) and supported by the National Institute for Health Research University College London Hospitals Biomedical Research Centre (to U. Menon and S.A. Gayther). USC: P01-CA17054, P30-CA14089, R01-CA61132, N01-PC67010, N01-CN025403 (to M.C. Pike), R03-CA113148, R03-CA115195 (to C.L. Pearce), and California Cancer Research Program (00-01389V-20170, 2II0200; to A. Wu); POL: Intramural Research Program of the NCI.