Genome-wide association analysis of more than 120,000 individuals identifies 15 new susceptibility loci for breast cancer.

Genome-wide association studies (GWAS) and large-scale replication studies have identified common variants in 79 loci associated with breast cancer, explaining ∼14% of the familial risk of the disease. To identify new susceptibility loci, we performed a meta-analysis of 11 GWAS, comprising 15,748 breast cancer cases and 18,084 controls together with 46,785 cases and 42,892 controls from 41 studies genotyped on a 211,155-marker custom array (iCOGS). Analyses were restricted to women of European ancestry. We generated genotypes for more than 11 million SNPs by imputation using the 1000 Genomes Project reference panel, and we identified 15 new loci associated with breast cancer at P < 5 × 10(-8). Combining association analysis with ChIP-seq chromatin binding data in mammary cell lines and ChIA-PET chromatin interaction data from ENCODE, we identified likely target genes in two regions: SETBP1 at 18q12.3 and RNF115 and PDZK1 at 1q21.1. One association appears to be driven by an amino acid substitution encoded in EXO1.

Breast cancer is the most common cancer in women worldwide[1]. The disease aggregates in families, and has an important inherited component. This inherited component is driven by a combination of rare variants, notably in BRCA1, BRCA2, PALB2, ATM and CHEK2 conferring a moderate or high lifetime risk of the disease, together with common variants at more than 70 loci, identified through GWAS and large scale replication studies[2-20]. Taken together, these loci explain approximately one-third of the excess familial risk of breast cancer.

The majority of susceptibility SNPs has been identified through the Breast Cancer Association Consortium (BCAC), a collaboration involving more than 50 case-control studies. We recently reported the results of a large-scale genotyping experiment within BCAC, which utilised a custom array (iCOGS) designed to study variants of interest for breast, ovarian and prostate cancers. iCOGS comprised more than 200,000 variants, of which 29,807 had been selected from combined analysis of nine breast cancer GWAS involving 10,052 breast cancer cases and 12,575 controls of European ancestry. In total, 45,290 breast cancer cases and 41,880 controls of European ancestry from 41 studies were genotyped with iCOGS, leading to the discovery of 41 novel susceptibility loci[16]. A parallel analysis identified four loci specific to oestrogen receptor (ER)-negative disease[17]. However, additional susceptibility loci may have been missed because they were not selected from the original GWAS, or not included on the array.

Genotype imputation is a powerful approach to infer missing genotypes using the genetic correlations defined in a densely genotyped reference panel, thus providing the opportunity to identify novel susceptibility variants even if not directly genotyped[21]. In this analysis we aimed to identify additional breast cancer susceptibility loci by utilising data from all 200k variants on the iCOGS array, and used imputation to estimate genotypes for more than 11M SNPs. We applied the same approach to data from 11 GWAS. After quality control (QC) exclusions, the dataset comprised 15,748 breast cancer cases and 18,084 controls from GWAS, and 46,785 cases and 42,892 controls from 41 studies genotyped with iCOGS (see Online Methods and Supplementary Tables 1a–1e). All subjects were women of European ancestry.

We imputed genotypes using the 1000 Genomes Project March 2012 release as the reference dataset (see Online Methods) The main analyses were based on ~11.6M SNPs that were imputed with imputation r2 >0.3 and had MAF>0.005 in at least one of the datasets[22].

Of common SNPs (MAF>0.05), 88% were imputed from the iCOGS array with r2>0.5; this compared to 99% of variants for the largest GWAS (UK2), which was genotyped using a 670k SNP array (Figure 1a and 1b, Supplementary Table 2). Thirty-seven per cent of common SNPs were imputed on the iCOGS with r2>0.9, compared with 85% for UK2. Thus, despite being designed as a follow-up of GWAS for different diseases rather than a genome-wide array, the majority of common variants could be imputed using the iCOGS, but the overall imputation quality was, poorer that from a standard GWAS array. Imputation quality decreased with decreasing allele frequency (Figure 1c and 1d, Supplementary Table 2).

Figure 1

Histograms of the imputation r2 a) Histogram of the imputation r2 for the iCOGS for variants with MAF>0.05 b) Histogram of the imputation r2 for the UK2 GWAS for variants with MAF>0.05 c) Histogram of the imputation r2 for the iCOGS for variants with MAF<=0.05 d) Histogram of the imputation r2 for the UK2 GWAS for variants with MAF<=0.05.

Log odds ratio estimates and standard errors were calculated for each dataset using logistic regression, adjusting for principal components where it was found to reduce substantially the inflation factor. We then combined the results from each dataset for variants with MAF >0.5% using a fixed effects meta-analysis[23]. More than 7,000 variants with a combined P<5×10−8 for association were identified, the large majority of which was in regions previously shown to be associated with breast cancer susceptibility. Of the 79 previously published breast cancer susceptibility loci identified in women of European ancestry, all but eight show evidence of association at P<5×10−8 for overall, ER-positive or ER-negative disease risk (Supplementary Tables 3a, 3b and 3c). For four of the eight variants, (rs1550623 on 2q31, rs11571833 on 13q13.1, rs12422552 on 12p13.1 and rs11242674 on 6p25.3), slightly weaker evidence of association was observed. One reported variant, rs7726159 did not reach P<5×10−8 in this (P=0.0017) or the previous analysis – it was identified through fine-mapping of the TERT region on 5p15.33[18]. One other variant in AKAP9, rs6964587 reported previously[19] did not reach P<5×10−8 but an alternative correlated with it did (P=3.67×10−8 for chr7:91681597:D; r2 between the two markers = 0.98). The two remaining variants (rs2380205 on 10p15 and rs1045485 at CASP8) were reported in earlier analysis[9,24] but did not even reach P<0.0001, suggesting that they may have been false positive reports. An alternative variant at CASP8, rs1830298 (r2=0.06, D’=1 with rs1045485 in 1000G CEU) did reach P<5×10−8 in this dataset[25].

To assess evidence for additional susceptibility loci, we removed all SNPs within 500kb of susceptibility variants identified previously in women of European ancestry[2-14,16-19], leaving 314 variants from 27 regions associated with breast cancer at P<5×10−8 (Supplementary Figures 1 and 2). The strongest associations were observed in a 610kb (b37 28,314,612- 28,928,858) interval on chromosome 22 (smallest P=8.2×10−22, for rs62237573). This interval lies approximately 100kb centromeric to CHEK2, and further analysis revealed that the associated SNPs were correlated with the CHEK2 founder variant 1100delC (strongest correlation r2=0.39 for SNP rs62235635), CHEK2 1100delC is known to be associated with breast cancer through candidate gene analysis, but has not previously generated an association in GWAS [26,27]. We performed an analysis adjusting for CHEK2 1100delC using data on ~40,000 samples that had been genotyped for this variant. The strongest associated variant in this subset was rs140914118; after adjustment for 1100delC the statistical significance diminished markedly (P=3.1×10−9 to P=0.78; Supplementary Figures 3a and 3b), suggesting that this signal is driven by CHEK2 1100delC.

Variants in four regions (DNAJC1, 5p12, PTHLH and MKL1) lay within 2Mb of a previously published susceptibility-associated SNP. In each case, these associations became weaker (no longer P<5×10−8) after adjustment for the previously associated SNP(s) in the region (data not shown). For four other regions, the significant variants were identified in just one GWAS, and failed imputation (r2<0.3) in the remaining datasets, including iCOGS; we did not consider these variants further.

To confirm the results for the remaining 18 regions, we performed re-imputation in the iCOGS dataset without phasing (See Online Methods). Fifteen loci remained associated with breast cancer at P<5×10−8 (Table 1 and Supplementary Table 4). For three of the loci, the most significant SNP, or a highly correlated SNP, had been directly genotyped on iCOGS (Supplementary Table 5); one, rs11205277, had been included on the array because it is associated with adult height[28], while the other two were selected based on evidence from the combined breast cancer GWAS but failed to reach genome-wide significance in the earlier analyses. We attempted to genotype the 12 remaining variants on a subset of ~4K samples to confirm the quality of the imputation (10 variants could be directly genotyped, for one region an alternative correlated variant was selected (Supplementary Table 5). For the 11 variants that could be assessed, the r2 between the observed and imputed genotypes were close to the r2 estimated in the imputation. Furthermore, the estimated effect sizes in the subset of individuals that we genotyped were similar to those obtained from the imputed genotypes (Supplementary Table 5). These results indicate that the analyses based on imputed genotype data were reliable.

Table 1

Results for the 15 regions with combined P<5×10−8. Results are shown for the strongest associated variant in the region.

Best variant	Locus	Position2	Alleles3	EAF4	r25	GWAS OR(95% CI)6	GWAS P7	iCOGS OR(95% CI)	iCOGS P	CombinedGWAS +iCOGS P	Genes within+/−2kb	Enhancers inMCF7/HMEC	eQTLs
rs12405132	1q21.1	145644984	C/T	0.36	0.96	0.96 (0.92–0.99)	0.00962	0.95 (0.93–0.97)	2.34×10−7	7.92×10−9	LOC10028814, NBPF10, RNF115	RNF115, POLR3C,PDZK1, PIAS3	-
rs12048493	1q21.2	149927034	A/C	0.34	0.76	1.04 (0.99–1.09)	0.121	1.07 (1.05–1.10)	1.66×10−9	1.10×10−9	-	-	-
rs72755295	1q43	242034263	A/G	0.03	0.94	1.19 (1.03–1.39)	0.021	1.15 (1.09–1.22)	2.60×10−7	1.82×10−8	EXO1	-	-
rs6796502	3p21.3	46866866	G/A	0.09	0.91	0.92 (0.87–0.98)	0.00657	0.92 (0.89–0.95)	8.13×10−7	1.84×10−8	-	-	-
rs13162653	5p15.1	16187528	G/T	0.45	0.72	0.92 (0.88–0.95)	5.18×10−6	0.95 (0.93–0.97)	1.71×10−6	1.08×10−10	-	-	-
rs2012709	5p13.3	32567732	C/T	0.46	0.81	1.06 (1.02–1.09)	0.00101	1.05 (1.03–1.08)	1.66×10−6	6.38×10−9	-	-	-
rs7707921	5q14	81538046	A/T	0.23	0.88	0.94 (0.9–0.98)	0.00302	0.93 (0.91–0.95)	4.09×10−9	5.00×10−11	ATG10	-	RPS23, ATP6AP1L
rs9257408	6p22.1	28926220	G/C	0.38	0.92	1.05 (1–1.1)	0.0372	1.05 (1.03–1.08)	4.53×10−7	4.84×10−8	-	-	-
rs4593472	7q32.3	130667121	C/T	0.35	1.00	0.92 (0.88–0.96)	2.57×10−5	0.95 (0.94–0.97)	3.97×10−6	1.83×10−9	FLJ43663	-	-
rs13365225	8p11.23	36858483	A/G	0.17	0.94	0.89 (0.85–0.93)	6.32×10−7	0.95 (0.93–0.98)	0.000159	1.06×10−8	-	-	-
rs13267382	8q23.3	117209548	G/A	0.36	0.97	1.07 (1.03–1.12)	0.000537	1.05 (1.03–1.07)	4.87×10−6	1.72×10−8	LINC00536	-	-
rs11627032	14q32.12	93104072	T/C	0.26	0.73	0.94 (0.9–0.98)	0.00114	0.94 (0.92–0.96)	1.06×10−6	4.48×10−9	RIN3	-	-
chr17:29230520	17q11.2	29230520	GGT/G	0.20	0.77	0.94 (0.89–0.98)	0.009	0.93 (0.91–0.96)	1.11×10−6	3.34×10−8	ATAD5	-	-
rs745570	17q25.3	77781725	A/G	0.50	0.93	0.94 (0.91–0.98)	0.000754	0.95 (0.93–0.97)	4.52×10−7	1.40×10−9	-	-	-
rs6507583	18q12.3	42399590	A/G	0.07	0.96	0.91 (0.85–0.98)	0.00803	0.91 (0.88–0.95)	1.21×10−6	3.20×10−8	SETBP1	SETBP1	-

Chromosome

Build 37 position

Reference/effect allele, based on the forward strand

Mean effect allele frequency over all controls

Imputation r2 in the iCOGS samples (calculated by the average info score from IMPUTEv2)

Per allele odds ratio for the minor allele relative to the major allele

P value for the 1df trend test

There was little or no evidence of heterogeneity in the per-allele odds ratios (ORs) among studies genotyped using iCOGS (Supplementary Table 6 and Supplementary Figure 4). There was little evidence for departure from a log-additive model for any locus, except for a borderline departure for rs6796502 (P=0.049) for which the ORs for heterozygotes and homozygotes for the risk associated allele were similar (Supplementary Table 6).

The estimated ORs for invasive versus in-situ disease were similar for all the loci (P>0.05) (Supplementary Table 7). For four of the loci, rs12405132, rs12048493, rs4593472 and rs6507583 the association was stronger for ER positive disease (case only P<0.05) (Supplementary Table 8). Seven of the loci were associated with ER-negative disease (P<0.05) but none had a stronger association for ER-negative than ER-positive disease. Two of the loci showed significant trends in the OR by age at diagnosis: for rs13162653, the OR was higher at younger ages (P=0.007), while for rs6507583, the OR was higher at older ages (P=0.006) (Supplementary Table 9). One of the variants, chr17:29230520:D in ATAD5 is correlated with a variant that has also been shown to be associated with serous ovarian cancer in a meta-analysis[29] (r2=0.93 between chr17:29230520:D and chr17:29181220:I).

To approach the task of identifying the likely causal variants and genes underlying these associations, we first defined the set of all SNPs correlated with each of the 15 lead SNPs and that could not be ruled out as potentially causal (based on a likelihood ratio 100:1[30]), resulting in a subset of 522 variants (Supplementary Table 10). One of the variants, rs72755295, lies in an intron of EXO1, encoding a protein involved in mismatch repair. It is strongly correlated with only one other variant, rs4149909, coding for an amino-acid substitution in EXO1 (p.Asn279Ser; CADD score 33[31]), suggesting that this variant is likely to be functionally related to breast cancer risk. None of the remaining SNPs lay within gene coding sequences, consistent with previous observations that most common cancer susceptibility variants are regulatory. For each of the remaining 520 variants, we then looked for enhancer elements in mammary cell lines, based on ENCODE ChIP-Seq data[32,33]. To identify potential gene targets, we combined this information with ENCODE ChIA-PET chromatin interaction data. We identified two regions in which the associated variants overlapped with putative enhancer sequences and for which consistent promoter interactions were predicted (Table 1). For rs12405132 at 1q21.1, we identified four potential interacting genes, RNF115, POLR3C, PDZK1 and PIAS3 (Figure 2). Of these, the strongest evidence was for RNF115 and PDZK1; three of the 64 potentially causal variants lay in interacting enhancer regions. RNF115 (also known as BCA2) is an E3 ubiquitin ligase RING finger protein that is overexpressed in ER-positive breast cancers[34]. PDZK1 is a scaffold protein that connects plasma membrane proteins and regulatory components, regulating their surface expression in epithelial cells apical domains, and has been proposed to act as an oncogene in breast cancer[35].

Figure 2

The chromosome 1 locus tagged by rs12405132 a) The Manhattan Plot displays the strength of genetic association (−log10 P) versus chromosomal position (Mb), where each dot presents a genotyped (solid black dot) or imputed (red circle) SNP (in the iCOGS stage). The purple horizontal line represents the threshold for genome-wide significance (P=5×10−8). Gene structures are depicted as well as the location of SNPs with MAF>0.01 which were neither imputed reliably nor genotyped. b) Mammary cell enhancer locations as defined in Corradin et al.[32], and Hnisz et al.[33], are shown where elements overlapping the best associated SNPs are labelled with their predicted target genes. A subset of ChiA-PET interactions in MCF7 cells (mediated by either RNApolII or ERa) between enhancers and their target gene promoters are also shown.

SNPs correlated with rs6507583 at 18q12.3 lay in regions interacting with the promoter of SETBP1 (Supplementary Figure 5). The encoded protein has been shown to bind the SET nuclear oncogene which is involved in DNA replication.

We utilised data from TCGA to assess associations between the 15 novel susceptibility variants and expression of neighbouring genes in breast tumors and normal breast tissue. One SNP, rs7707921, was strongly associated with RPS23 expression in all tissues (Supplementary Table 11, Supplementary Figure 6). However, stronger associations with expression were observed with more telomeric SNPs that were less strongly associated with disease risk (top eQTL SNP rs3739: P=10−23, P-risk=5.28×10−7), suggesting that this association may be coincidental. SNP, rs7707921 was also more weakly associated with expression of ATP6AP1L (P=5.6×10−5 in tumours, P=0.066 in normal tissue).

Based on the estimated ORs in the iCOGS stage (all but one of which were in the range 1.05–1.10), the 15 novel loci identified here would explain a further ~2% of the 2-fold familial risk of breast cancer. Taken together with previously identified loci, more than 90 independent common susceptibility loci for breast cancer have been identified, explaining ~16% of the familial risk. We estimate assuming a log-additive model that, based on genotypes for variants at these loci, approximately 5% of women in the general population have a >2 fold increased risk of breast cancer and 0.7% of women have a >3 fold increased risk. In the current analyses, more than 50% of variants with MAF>0.005 in subjects of European ancestry were well imputable (r2>0.5) These results suggest that, while there may be further susceptibility variants with comparable associated effects that were not well imputed, the identification of many additional loci will require larger association studies. In the meantime, inclusion of these additional loci in polygenic risk scores will improve our ability to discriminate between high and low risk individuals, potentially improving breast cancer screening and prevention.

Online Methods

Details of the subjects, genotyping and QC measures for the GWAS and iCOGS data are described elsewhere[12,14,16,36,37]. All participating studies were approved by their appropriate ethics review board and all subjects provided informed consent. Analyses were restricted to women of European ancestry. All imputations were performed using the 1000 Genomes Project March 2012 release as the reference panel. Of the 11 GWAS, 8 (C-BCAC) plus a subset of the BPC3 GWAS (CGEMS) were used in the combined GWAS analysis that nominated 29,807 SNPs for the array. The BPC3 and TNBCC GWAS nominated additional SNPs with evidence for association with ER-negative or triple-negative (ER-, PR- and HER2- negative) breast cancer. The EBCG GWAS was not used to nominate SNPs for the iCOGS array.

For eight GWAS (C-BCAC), genotypes were imputed in a two-stage procedure, using SHAPEIT to derive phased genotypes and IMPUTEv2 to perform the imputation on the phased data [22]. We performed the imputation using 5Mb non-overlapping intervals for the whole genome. OR estimates and standard errors where obtained using logistic regression with SNPTEST [21]. For two of the studies we adjusted for the 3 leading principal components as it was found to reduce materially the inflation factor; for the rest of the studies no such adjustment was necessary. For the remaining three GWAS (BPC3, TNBCC and EBCG), imputation was performed using MACH and Minimac[23]. Genomic control adjustment was applied to each GWAS as previously described[16]. The iCOGS data were also imputed in a two-stage procedure using SHAPEIT and IMPUTEv2, again using 5Mb non-overlapping intervals. We split the ~90K samples into 10 subsets, where possible keeping subjects from the same study in the same subset. We obtained OR estimates and standard errors using logistic regression adjusting for study and 9 principal components.

For the regions showing evidence of association we repeated the imputation in iCOGS, using IMPUTEv2 but without pre-phasing in SHAPEIT to improve imputation accuracy. We also increased the number of MCMC iterations from 30 to 90, and increased the buffer region from 250kb to 500kb.

Meta-analysis

OR estimates and standard errors were combined in a fixed effects inverse variance meta-analysis using METAL[23]. For the GWAS, results were included in the analysis for all SNPs with MAF>0.01 and imputation r2>0.3, except for the TN GWAS where the criteria were r2>0.9 and MAF>0.05. For iCOGS, we included all SNPs with r2>=0.3 and MAF>0.005.

Confirmatory genotyping

The best variant in each region after the re-imputation and meta-analysis was genotyped in 4123 samples from SEARCH, using Taqman according to the manufacturer’s instructions. The squared correlations between the observed genotypes and the genotypes estimated by imputation are shown in Supplementary Table 5. For all the imputed SNPs the squared correlations was greater than 0.7, the call-rates were >=0.98 and there was no evidence of departure of genotype frequencies from those expected under HWE (p>0.1).

eQTL analyses

Germline genotype, mRNA expression, and somatic copy number data for samples taken from breast tumours and tumour-adjacent normal tissue were obtained from The Cancer Genome Atlas[38]. The copy number and genotype data were measured using the Affymetrix Genome-Wide Human SNP 6.0 platform. For the mRNA expression data, we used the expression profiles obtained using the Agilent G4502A-07-3 microarray. The genotype data were subjected to the following quality control filters. SNPs were excluded in case of low frequency (MAF < 1%), low call rate (< 95%,) or departure from Hardy-Weinberg equilibrium at P < 1 × 1013. Individuals were excluded based on low call rate (< 95%), or high heterozygosity (false discovery rate < 1%). Furthermore, individuals were also excluded in case of non-European ancestry, or male gender. Quality control and intersection with the other genomic data types resulted in 380 tumour samples and 56 normal samples.

The genotype data were imputed as described above. eQTL analysis was performed using linear regression with SNPTEST, regressing the mRNA expression of selected candidate genes on the imputed genotype. For each gene, we performed the eQTL analysis against every microarray probe that uniquely maps to that gene. We adjusted the analyses for somatic copy number of the gene, and for SNPs that intersect the probe sequence, provided that their MAF exceeds 1% in individuals of European ancestry in the 1,000 Genomes data.

Enhancer analyses

Maps of enhancer regions with predicted target genes were obtained from Hnisz et al.[33], and Corradin et al.[32]. Enhancers active in the mammary cell types MCF7, HMEC and HCC1954 were intersected with candidate causal variants using Galaxy. ENCODE ChIA-PET chromatin interaction data from MCF7 cells (mediated by RNApolII and ERα) were downloaded using the UCSC Table browser. Galaxy was used to identify ChIA-PET interactions between an implicated mammary cell enhancer (containing a strongly associated variant) and a predicted gene promoter (defined as regions 3 kb upstream and 1 kb downstream of the transcription start site).