Identification of six new susceptibility loci for invasive epithelial ovarian cancer.

We obtained data on BRCA1 and BRCA2 mutation carriers through CIMBA. Eligibility in CIMBA is restricted to females 18 years or older with pathogenic mutations in BRCA1 or BRCA2 . The majority of the participants were sampled through cancer genetics clinics 15 , including some related participants. Fifty-four studies from 27 countries contributed data. After quality control, data were available on 15,252 BRCA1 mutation carriers and 8,211 BRCA2 mutation carriers, of whom 2,462 and 631, respectively, were affected with EOC ( Supplementary Table 1 ). Data were available for the stage 1 of three population-based EOC GWAS. These included 2,165 cases and 2,564 controls from a GWAS from North America (“US GWAS”) 39 , 1,762 cases and 6,118 controls from a UK-based GWAS (“UK GWAS”) 6 , and 441 cases and 441 controls from the Mayo GWAS. Furthermore, 11,069 cases and 21,722 controls were genotyped using the iCOGS array (“OCAC-iCOGS” stage data). Overall, 43 studies from 11 countries provided data on 15,347 women diagnosed with invasive epithelial EOC, 9,627 of whom were diagnosed with serous EOC, and 30,845 controls from the general population. All subjects included in this analysis were of European descent and provided written informed consent as well as data and blood samples under ethically approved protocols. Further details of the OCAC and CIMBA study populations as well as the genotyping, quality control and statistical analyses have been described elsewhere 7 , 13 , 16 . Genotyping and imputation details for each study are shown in Supplementary Table 1 . To evaluate the accuracy of the imputation of the SNPs we found to be associated with EOC risk, we genotyped rs17329882 (4q26) and rs635634 (9q34.2) in a subset of 3,541 subjects from CIMBA using Sequenon’s iPLEX technology. The lead SNP at 17q11.2, chr17:29181220:I failed iPLEX design. We performed quality control of the iPLEX data according to the CIMBA guidelines. After quality control, we used the imputation results to generate the expected allele dosage for each genotyped sample and computed the Pearson product-moment correlation coefficient between the expected allele dosage and the observed genotype. The squared correlation coefficient was compared to the imputation accuracy as estimated from the imputation. We carried out quality control separately for BRCA1 carriers, BRCA2 carriers, the three OCAC GWAS, and OCAC-iCOGS samples, but quality criteria were mostly consistent across studies. We excluded samples if they were not of European ancestry, if they had a genotyping call rate < 95%, low or high heterozygosity, if they were not female or had ambiguous sex, or were duplicates (cryptic or intended). In OCAC studies, one individual was excluded from each pair of samples found to be first-degree relatives and duplicate samples between the iCOGS stage and any of the GWAS were excluded from the iCOGS data. SNPs were excluded if they were monomorphic, had call rate<95%, showed evidence of deviation from Hardy-Weinberg equilibrium or had low concordance between duplicate pairs. For the Mayo GWAS and the UK GWAS, we also excluded rare SNPs (MAF<1% or allele count <5, respectively). We visually inspected genotype cluster plots for all SNPs with P<10 −5 from each of the newly identified loci. We used the R GenABEL library version 1.6.7 for quality control 40 . Genotype data were available for analysis from iCOGS for 199,526 SNPs in OCAC-iCOGS, 200,720 SNPs in BRCA1 mutation carriers, and 200,908 SNPs in BRCA2 mutation carriers. After QC, for the GWAS, data were available on 492,956 SNPs for the US GWAS, 543,529 SNPs for the UK GWAS and 1,587,051 SNPs for the Mayo GWAS ( Supplementary Table 2 ). We performed imputation separately for BRCA1 carriers, BRCA2 carriers, OCAC-iCOGS samples and each of the OCAC GWAS. We imputed variants from the 1000 Genomes Project data using the v3 April 2012 release 17 as the reference panel. For OCAC-iCOGS, the UK GWAS and the Mayo GWAS, imputation was based on the 1000 Genomes Project data with singleton sites removed. To improve computation efficiency we initially used a two-step procedure, which involved pre-phasing in the first step and imputation of the phased data in the second. We carried out pre-phasing using the SHAPEIT software 41 . We used the IMPUTE version 2 software for the subsequent imputation 42 for all studies with the exception of the US GWAS for which the MACH algorithm implemented in the minimac software version 2012.8.15, mach version 1.0.18 was used. To perform the imputation we divided the data into segments of approximately 5Mb each. We excluded SNPs from the association analysis if their imputation accuracy was r 2 <0.3 or their minor allele frequency (MAF) was <0.005 in BRCA1 or BRCA2 carriers or if their accuracy was r 2 <0.25 in OCAC-iCOGS, the UK GWAS, UK GWAS or Mayo GWAS. We performed more accurate imputation for the regions around the novel EOC loci from the joint analysis of the data from BRCA1 and BRCA2 carriers and the general population (any SNP with P<5×10 −8 ). The boundaries of these regions were set +/− 500kb from any significantly associated SNP in the region. As in the first run, the 1000 Genomes Project data v3 were used as the reference panel and the software IMPUTE2 was applied. However, for the second round of imputation, we imputed genotypes without pre-phasing in order to improve accuracy. To further increase the imputation accuracy we changed some of the default parameters in the imputation procedure. These included an increase of the MCMC iterations to 90 (out of which the first 15 were used as burn-in), an increase of the buffer region to 500kb and an increase of the number of haplotypes used as templates when phasing observed genotypes to 100. These changes were applied consistently for all data sets. We evaluated the association between genotype and disease using logistic regression by estimating the associations with each additional copy of the minor allele (log-additive models). The analysis was adjusted for study and for population substructure by including the eigenvectors of the first five ancestry specific principal components as covariates in the model. We used the same approach to evaluate the SNP associations with serous ovarian cancer after excluding all cases with any other or with unknown tumour subtype. For imputed SNPs we used expected dosages in the logistic regression model to estimate SNP effect sizes and p-values. We carried out analyses separately for OCAC-iCOGS and the three GWAS and pooled thereafter using a fixed effects meta-analysis. We carried out the analysis of re-imputed genotypes of putative novel susceptibility loci jointly for the OCAC-iCOGS and GWAS samples. All results are based on the combined data from iCOGS and the three GWAS. We used custom written software for the analysis. We carried out the ovarian cancer association analyses separately for BRCA1 and BRCA2 mutation carriers. The primary analysis was carried out within a survival analysis framework with time to ovarian cancer diagnosis as the endpoint. Mutation carriers were followed until the age of ovarian cancer diagnosis, or risk-reducing salpingo-oophorectomy (RRSO) or age at last observation. Breast cancer diagnosis was not considered as a censoring event. In order to account for the non-random sampling of BRCA1 and BRCA2 mutation carriers with respect to their disease status we conducted the analyses by modelling the retrospective likelihood of the observed genotypes conditional on the disease phenotype 18 . We assessed the associations between genotype and risk of ovarian cancer using the 1 degree of freedom score test statistic based on the retrospective likelihood 18 , 43 . To account for the non-independence among related individuals in the sample, we used an adjusted version of the score test statistic, which uses a kinship adjusted variance of the score 44 . We evaluated associations between imputed genotypes and ovarian cancer risk using a version of the score test as described above but with the posterior genotype probabilities replacing the genotypes. All analyses were stratified by the country of origin of the samples. We carried out the retrospective likelihood analyses in CIMBA using custom written functions in Fortran and Python. The score test statistic was implemented in R version 3.0.1 45 . We evaluated whether there is evidence for multiple independent association signals in the region around each newly identified locus by evaluating the associations of genetic variants in the region while adjusting for the SNP with the smallest meta-analysis p-value in the respective region. This was done separately for BRCA1 carriers, BRCA2 carriers and OCAC. For one of the novel associations, it was not possible to confirm the imputation accuracy of the lead SNP chr17:29181220:I at 17q11.2 through genotyping. Therefore, we inferred two-allele haplotypes for rs9910051 and rs3764419, highly correlated with the lead SNP (r 2 =0.95), using an in-house program. These variants were genotyped on the iCOGS array and therefore this analysis was restricted to 14,733 ovarian cancer cases and 9,165 controls from OCAC-COGS, and 8,185 BRCA2 mutation carriers that had available genotypes for both variants based on iCOGS. The association between the AA haplotype and risk was tested using logistic regression in OCAC and using Cox regression in BRCA2 mutation carriers. We conducted a meta-analysis of the EOC associations in BRCA1, BRCA2 carriers and the general population for genotyped and imputed SNPs using an inverse variance approach assuming fixed effects. We combined the logarithm of the per-allele hazard ratio estimate for the association with EOC risk in BRCA1 and BRCA2 mutation carriers and the logarithm of the per-allele odds ratio estimate for the association with disease status in OCAC. For the associations in BRCA1 and BRCA2 carriers, we used the kinship adjusted variance estimator 44 which allows for inclusion of related individuals in the analysis. We only used SNPs with results in OCAC and in at least one of the BRCA1 or the BRCA2 analyses. We carried out two separate meta-analyses, one for the associations with EOC in BRCA1 carriers, BRCA2 carriers and EOC in OCAC, irrespective of tumour histological subtype, and a second using only the associations with serous EOC in OCAC. The number of BRCA1 and BRCA2 samples with tumour histology information was too small to allow for subgroup analyses. However, previous studies have demonstrated that the majority of EOCs in BRCA1 and BRCA2 mutation carriers are high-grade serous 49 – 53 . Meta-analyses were carried out using the software “metal”, 2011-03-25 release 54 . In order to identify a set of potentially causal variants we excluded SNPs with a likelihood of being causal of less than 1:100, by comparing the likelihood of each SNP from the association analysis with the one of the most strongly associated SNP 46 . The remaining variants were then analysed using pupasuite 3.1 to identify potentially functional variants ( Supplementary Table 9 ). Early-passage primary normal ovarian surface epithelial cells (OSECs) and fallopian tube epithelial cells were harvested from disease-free ovaries and fallopian tubes. Normal ovarian epithelial cells were collected by brushing the surface of the ovary with a sterile cytobrush, and were cultured in NOSE-CM 55 . Fallopian tube epithelial cells were harvested by Pronase digestion as previously described 56 , plated onto collagen-coated plastics (Sigma) and cultured in DMEM/F12 (Sigma-Aldrich) supplemented with 2% Ultroser G (BioSepra) and 1× penicillin/streptomycin (Lonza). By the time of RNA harvesting, fallopian tube cultures tested consisted of PAX8 positive fallopian tube secretory epithelial cells (FTSECs), consistent with previous observations that ciliated epithelial cells from the fallopian tube do not proliferate in vitro . For gene expression analysis, RNA was harvested from 59 early passage samples: 54 OSECs and 5 FTSECs from cell cultures harvested at ~80% confluency using the QIAgen miRNAeasy kit with on-column DNase 1 digestion. 500ng RNA was reverse transcribed using the Superscript III kit (Life Technologies). We preamplified 10ng cDNA using the TaqMan ® Preamp Mastermix; the resulting product was diluted 1:60 and used to quantify gene expression using the following TaqMan ® gene expression probes: WNT4, Hs01573504_m1; RSPO1, Hs00543475_m1; SYNPO2, Hs00326493_m1; ATAD5, Hs00227495_m1 and GPX6, Hs00699698_m1. Four control genes were also included: ACTB, Hs00357333_g1; GAPDH, Hs02758991_g1; HMBS, Hs00609293_g1 and HPRT1 Hs02800695_m1 (all Life Technologies). Assays were run on an ABI 7900HT Fast Real-Time PCR system (Life Technologies). Expression levels for each gene were normalized to the average of all four control genes. Relative expression levels were calculated using the δδCt method. Genotyping was performed on the iCOGs chips, as described above. Where genotyping data were not available for the most risk-associated SNP, the next most significant SNP was used: rs3820282 at 1p36, rs12023270 at 1p34.3, rs752097 at 4q26, rs445870 at 6p22.1, rs505922 at 9q34.2 and rs3764419 at 17q11.2. Correlations between genotype and gene expression were calculated in ‘R’. Genotype specific gene expression in the normal tissue cell lines (eQTL analysis) was compared using the Jonckheere-Terpstra test. IData were normalized to the four control genes and we tested for eQTL associations, grouping OSECs and FTSECs together. Secondly, OSECs were analysed alone. eQTL analyses were performed using 3 genotype groups, or two groups (with the rare homozygote samples grouped together with the heterozygote samples). eQTL analysis in primary tumours was based on the publicly available data available from The Cancer Genome Atlas (TCGA) project, which includes 489 primary high grade serous ovarian cancers. The methods have been described elsewhere 57 . Briefly, we determined the ancestry for each case based on the germ line genotype data using EIGENSTRAT software with 415 HapMap genotype profiles as a control set. Only populations of Northern and Western European ancestries were included. We first performed a cis -eQTL analyses using a method we described previously, in which the association between 906,600 germline genotypes and the expression levels of mRNA or miRNA (located within 500Kb on either side of the variant) were evaluated using linear regression model with the effects of somatic copy number and CpG methylation being deducted (For miRNA expression, the effect of CpG methylation is not adjusted for since the data are not available). To adjust for multiple tests, we adjusted the test P values using Benjamini-Hochberg method. A significant association was defined by a false discovery rate (FDR) of less than 0.1. Having established a genome-wide cis -eQTL associaitions in this series of tumours, we then evaluated cis -eQTL associations for the top risk associations between each of the six new loci and the gene in closest proximity to the risk SNP. For each risk locus, we retrieved the genotype of all SNPs in ovarian cancer cases based on the Affymetrix 6.0 array. Using these genotypes and the impute2 March 2012 1000 Genomes Phase I integrated variant cosmopolitan reference panel of 1,092 individuals (Haplotypes were phased via SHAPEIT), we imputed the genotypes of SNPs in the 1000 Genomes Project in the target regions for TCGA samples 58 . For each risk locus where data for the most risk-associated variant were not available, we retrieved the imputed variants tightly correlated with the most risk-associated variant. We then tested for association between imputed SNPs and gene expression using the linear regression algorithm described above, where each imputed SNP was coded as an expected allele count. Again, significant associations are defined by a false discovery rate (FDR) of less than 0.1. We performed genome-wide formaldehyde assisted regulatory element (FAIRE) and ChIP seq with histone 3 lysine 27 acetylation (H3K27ac) and histone 3 lysine 4 monomethylation (H3K4me) for two normal OSECs, two normal FTSECs and two HGSOC cell lines (UWB1.289 and CAOV3) [Shen et al. in preparation]. These datasets annotate epigenetic signatures of open chromatin, and collectively indicate transcriptional enhancer regions. We analysed the FAIRE-seq and ChIP-seq datasets and publically available genomic data on promoter and UTR domains, intron/exon boundaries, and positions of non-coding RNA transcripts to identify SNPs from the 100:1 likely causal set that align with biofeatures that may provide evidence of SNP functionality. The Cancer Genome Atlas (TCGA) Project and COSMIC Datasets TCGA has performed extensive genomic analysis of tumours from a large number of tissue types including almost 500 high-grade serous ovarian tumours. These data include somatic mutations, DNA copy number, mRNA and miRNA expression and DNA methylation. COSMIC is the catalogue of somatic mutations in cancer that collates information on mutations in tumours from the published literature 59 . They have also identified The Cancer Gene Census, which is a list of genes known to be involved in cancer. Data are available on a large number of tissue types, including 2,809 epithelial ovarian tumours. We analysed all genes for coding somatic sequence mutations generated from either whole exome or whole genome sequencing. In TCGA, whole exome sequencing data were available for 316 high-grade serous EOC cases. In addition, we determined whether mutations had been reported in COSMIC 59 and whether the gene was a known cancer gene in the Sanger Cancer Gene Census. Normalized and gene expression values (Level 3) gene expression profiling data were obtain from the TCGA data portal for three different platforms (Agilent, Affymetrix HuEx and Affymetrix U133A). We analysed only the 489 primary serous ovarian tumour samples included in the final clustering analysis 58 and eight normal fallopian tube samples. The boxplot function in R was used to compare ovarian tumour samples to the fallopian tube for 91 coding genes with expression data on any platform within a 1MB region around the most significant SNP at the six loci. A difference in relative expression between EOC and normal tissue was carried out using the Wilcoxon rank-sum test. Serous EOC samples for 481 tumours with log2 copy number data were analysed using the cBio portal for analysis of TCGA data 60 , 61 . For each gene in a region the classes of copy number; homozygous deletion, heterozygous loss, diploid, gain, and amplification were queried individually using the advanced onco query language (OQL) option. The frequency of gain and amplification were combined as “gain”, and homozygous deletion and heterozygous loss were combined as “loss”. Serous EOC samples for 316 complete tumours (those with CNA, mRNA and sequencing data) were analysed. Graphs were generated using the cBio portal for analysis of TCGA data and the setting were mRNA expression data Z-score (all genes) with the Z-score threshold of 2 (default setting) and putative copy number alterations (GISTIC). The Z-score is the number of standard deviations away from the mean of expression in the reference population. GISTIC is an algorithm that attempts to identify significantly altered regions of amplification or deletion across sets of patients. The putative causal SNPs at the 1p36 locus lie in the WNT4 promoter and so we tested their effect on transcription in a luciferase reporter assay ( Fig. 2D ). Wild-type and risk haplotype (comprising five correlated variants) sequences corresponding to the region bound by hg19 co-ordinates chr1:22469416-22470869 were generated by Custom Gene Synthesis (GenScript Corporation), and then sub-cloned into pGL3-basic (Promega). Equimolar amounts of luciferase constructs (800 ng) and pRL-TK Renilla (50 ng) were co-transfected into ~8 × 10 4 iOSE4 62 normal ovarian cells in triplicate wells of 24 well plates using LipoFectamine 2000 (Life Technologies). Independent transfections were repeated three times. The Dual-Glo Luciferase Assay kit (Promega) was used to assay luciferase activity 24 hours post transfection using a BioTek Synergy H4 plate reader. The iOSE-4 cell line (derived by K. Lawrenson) was maintained under standard conditions and routinely tested for Mycoplasma and short tandem repeat profiled.