Genome-wide association analysis in East Asians identifies breast cancer susceptibility loci at 1q32.1, 5q14.3 and 15q26.1.

In a three-stage genome-wide association study among East Asian women including 22,780 cases and 24,181 controls, we identified 3 genetic loci newly associated with breast cancer risk, including rs4951011 at 1q32.1 (in intron 2 of the ZC3H11A gene; P=8.82×10(-9)), rs10474352 at 5q14.3 (near the ARRDC3 gene; P=1.67×10(-9)) and rs2290203 at 15q26.1 (in intron 14 of the PRC1 gene; P=4.25×10(-8)). We replicated these associations in 16,003 cases and 41,335 controls of European ancestry (P=0.030, 0.004 and 0.010, respectively). Data from the ENCODE Project suggest that variants rs4951011 and rs10474352 might be located in an enhancer region and transcription factor binding sites, respectively. This study provides additional insights into the genetics and biology of breast cancer.

Breast cancer is one of the most common malignancies among women worldwide. Genetic factors play a significant role in breast cancer etiology[1, 2]. To date, genome-wide association studies (GWAS) have identified approximately 75 genetic loci to be associated with breast cancer risk[2-5]. With the exception of the studies we have conducted among East Asian women[6-9] and one study conducted among women of African ancestry[10], all other published GWAS have been conducted among women of European ancestry. Genetic risk variants identified to date from GWAS explain only about 10% of familial risk for breast cancer in East Asian women[3]. Given the difference in genetic architecture and environmental exposures between women of European and East Asian ancestry, additional GWAS need to be conducted among East Asian women to fully uncover the genetic basis of breast cancer risk.

The current study was conducted as part of the Asia Breast Cancer Consortium (ABCC) to search for additional novel susceptibility loci for breast cancer. Included in this study are data obtained from 22,780 breast cancer cases and 24,181 controls who were recruited in 14 studies conducted in multiple Asian countries (Supplementary Table 1). The discovery stage (Stage 1) included two GWAS, in which 5,285 Chinese women (SBCGS-1) and 4,777 Korean women (SeBCS1) were scanned primarily using the Affymetrix Genome-Wide Human SNP Array 6.0, which consists of 906,602 single nucleotide polymorphisms (SNPs). After applying quality control (QC) filters described previously[6, 9, 11], 5,152 Chinese women (2,867 cases and 2,285 controls; 677,157 SNPs) and 4,298 Korean women (2,246 cases and 2,052 controls; 555,117 SNPs) remained in the current analysis. Imputation was conducted for each study following the MACH algorithm[12] using HapMap II release 22 CHB and JPT data (2,416,663 SNPs) as the reference. Only SNPs with a high imputation quality score (RSQR ≥ 0.50) were analyzed for associations with breast cancer risk. In the analyses of data from Chinese and Korean women, a total of 1,930,412 SNPs and 1,907,146 SNPs, respectively, were included. A meta-analysis of these GWAS data was conducted using a fixed-effects, inverse-variance meta-analysis with the METAL program[13]. There was little evidence for inflation in the association test statistics for the studies included in Stage 1 (genomic inflation factors (λ): λ= 1.0426 for SBCGS-1, λ= 1.0431 for SeBCS1, and λ = 1.0499 for both studies combined; Supplementary Figure 1). When scaled to a study of 1,000 cases and 1,000 controls, λ1,000 were 1.02, 1.02, and 1.01, respectively.

To select SNPs for the Stage 2 replication, we used the following criteria: (1) P < 0.05 in the Stage 1 meta results; (2) same direction of association in both Stage 1 studies; (3) no heterogeneity observed between the two Stage 1 studies (P > 0.05 and I < 25%); (4) an imputation score of RSQR > 0.5 in both Stage 1 studies; (5) a minor allele frequency (MAF) of > 0.05 in both Stage 1 studies; and (6) not in strong LD (r2 < 0.5) with any of the known breast cancer susceptibility loci or any SNPs we had evaluated previously[3, 6–9]. For SNPs that met the above criteria but were in LD (r2 > 0.5) with each other, we selected only one SNP for replication. A total of 4,598 SNPs were selected, and assays for 4,071 SNPs were successfully designed using Illumina Infinium assays as a part of large-scale genotyping effort. Of the 4,071 SNPs, 3,850 SNPs were successfully genotyped in an independent set of 3,944 cases and 3,980 controls selected from the Shanghai studies (SBCGS-2). After QC exclusions, 3,678 SNPs were included in the analyses of 3,472 cases and 3,595 controls.

For Stage 3, the top 50 SNPs were selected for further replication in an independent set of 14,195 cases and 16,249 controls from 10 studies participating in the ABCC based on the following criteria: (1) P < 0.005 in the meta-analysis of Stage 1 and 2 data and (2) same direction of association in both Stages 1 and 2. Of the 50 SNPs evaluated in Stage 3, 11 SNPs showed an association with breast cancer risk at P < 0.05 (Supplementary Table 2). Combined analyses of data from all three stages identified three SNPs that were associated with breast cancer risk at the genome-wide significance level (P < 5.0 × 10−8): rs4951011 at 1q32.1, odds ratio (OR) = 1.09, P = 8.82 × 10−9; rs10474352 at 5q14.3, OR = 1.09, P = 1.67 × 10−9; and rs2290203 at 15q26.1, OR = 1.08, P = 4.25 × 10−8 (Table 1). The association between breast cancer risk and each of these three SNPs was consistent across the studies included in the ABCC (Figure 1), and none of the tests for heterogeneity were statistically significant (P > 0.05) (Table 1). No significant heterogeneity was found for the association of these three SNPs with breast cancer risk among Chinese, Japanese, or Korean women (Supplementary Table 3). One additional SNP showed an association with breast cancer risk with a P-value near the conventional GWAS significance level (rs11082321 at 18q11.2, OR = 1.08, P = 6.77× 10−7) (Supplementary Table 2).

Table 1

Associations of breast cancer risk with newly identified risk variants: Results from the Asia Breast Cancer Consortium

SNP (Allelesa)	Frequencyb	Locus (Positionc)	Closest gene (Annotation)	Stage	Per-allele association		P for heterogeneityf
					OR (95% CI)d	Pe
rs4951011 (G/A)	0.282	1q32.1 (202,032,954)	ZC3H11A (Intron 2)	Stage 1	1.09 (1.02–1.17)	0.007
				Stage 2	1.10 (1.02–1.18)	0.011
				Stage 3	1.08 (1.05–1.12)	1.02 × 10−5
				Combined	1.09 (1.06–1.12)	8.82 × 10−9	0.98
rs10474352 (C/T)	0.482	5q14.3 (90,767,981)	ARRDC3 (Intergenic)	Stage 1	1.09 (1.03–1.17)	0.006
				Stage 2	1.12 (1.05–1.20)	7.06 × 10−4
				Stage 3	1.08 (1.04–1.12)	1.92 × 10−5
				Combined	1.09 (1.06–1.12)	1.67 × 10−9	0.50
rs2290203 (G/A)	0.504	15q26.1 (89,313,071)	PRC1 (Intron 14)	Stage 1	1.08 (1.02–1.14)	0.012
				Stage 2	1.19 (1.10–1.30)	4.97 × 10−5
				Stage 3	1.06 (1.03–1.10)	2.45 × 10−4
				Combined	1.08 (1.05–1.11)	4.25 × 10−8	0.06

Abbreviations: OR, odds ratio; CI, confidence interval.

Risk/reference allele; risk allele is shown in bold.

Risk allele frequency in controls from all three stages combined.

Chromosome position (bp) based on NCBI Human Genome Build 36.

Per-allele OR (95% CI) was adjusted for age and principal components in each study; summary OR (95% CI) was obtained using fixed-effect meta-analysis in each stage.

Derived from a weighted z-statistic–based meta-analysis.

P for heterogeneity across studies in all stages was calculated using Cochran’s Q test.

Figure 1

Forest plots for risk variants in the three newly identified breast cancer risk loci by study site and stage

Per-allele ORs are presented. The size of the box is proportional to the number of cases and controls in each study. (a): rs4951011, (b): rs10474352, (c): rs2290203.

The associations of SNPs rs10474352 and rs2290203 appeared to be stronger for estrogen receptor (ER)-positive breast cancer than for ER-negative breast cancer, and the heterogeneity test was of borderline significance for rs10474352 (P = 0.085) (Supplementary Table 4). The associations of rs4951011 with breast cancer risk were similar for ER-positive and ER-negative breast cancer.

We evaluated the three newly identified risk variants for associations with breast cancer risk in European-ancestry women using data from 16,003 cases and 41,335 controls derived from twelve breast cancer GWAS and included in the DRIVE GAME-ON Consortium. SNPs rs4951011, rs10474352, and rs2290203 were all associated with breast cancer risk in women of European ancestry at P < 0.05 with the same direction of association as observed in East Asian women (Supplementary Table 5). However, the strength of the associations was weaker in European-ancestry women than in East Asian-ancestry women, and the frequencies of the risk alleles were quite different between these two populations.

We evaluated and annotated putative functional variants and candidate genes in each of the three newly identified loci using data from the Encyclopedia of DNA Elements (ENCODE)[14], The Cancer Genome Atlas (TCGA) breast cancer project[15], expression quantitative trait locus (eQTL) databases[16] as well as RegulomeDB[17] and HaploReg v2[18]. We summarize the results below for each locus.

SNP rs10474352 is located on 5q14.3, 53,078 bp upstream of the ARRDC3 gene (Figure 2b). The ARRDC3 gene is a member of the arrestin gene family and is suspected of playing a role in breast cancer development. A gene cluster at 5q11-q23 that includes ARRDC3 was found to be deleted in 17% of breast cancer tumor tissue[19]. Up regulation of the ARRDC3 gene in a breast cancer cell line has been shown to repress cell proliferation, migration, invasion, and in vivo tumorigenesis[20]. We evaluated ARRDC3 gene expression in 87 breast cancer cases included in TCGA. The expression level of the ARRDC3 gene was significantly lower in tumor tissue than in adjacent normal tissue (P = 1.88 × 10−18) (Supplementary Table 6). This is consistent with a previous study showing that expression levels of the ARRDC3 gene were lower in breast tumor tissue compared with normal tissue and in metastatic tumor tissue compared with primary tumor tissue[20]. Furthermore, lower ARRDC3 expression in tumor tissue has been associated with poorer disease-free survival in breast cancer patients[20]. A search of RegulomeDB[17] and HaploReg[18] indicated that rs10474352 may be located in predicted AP-1and VDR motifs (Supplementary Table 7), suggesting a potential regulatory role. We evaluated whether SNPs in this locus are cis-eQTL for other genes by analyzing TCGA breast cancer data. Our analysis revealed no evidence that this SNP or its correlated SNPs are cis-eQTLs for any genes in this locus. Recently, a SNP located ~596 kb upstream of the ARRDC3 gene, rs421379, was found to be associated with prognosis for early-onset breast cancer in a GWAS[21]. However, rs421379 is not in LD with rs10474352 (r=0 in both ASN and CEU data), the SNP in close proximity to ARRDC3 that was identified in our study. Furthermore, in our study, rs421379 had a low MAF (0.03–0.04) and was not associated with breast cancer risk (P = 0.2484 in Stage 1).

Figure 2

Regional plots of association results for the three newly identified risk loci for breast cancer

For each plot, the −log10 (P values) (y axis) of SNPs are shown according to their chromosomal positions (x axis) in NCBI Build 36. The color of SNPs represents their LD (r, HapMap Asian) with the index SNP at each locus. With the exception of the index SNPs, which are shown as purple diamonds for Stage 1 and purple circles for the meta-analyses of all studies, data shown for all other SNPs are from Stage 1 only. (a): rs4951011 (1q32.1), (b): rs10474352 (5q14.3), (c): rs2290203 (15q26.1).

SNP rs2290203 is located in intron 14 of the protein regulator of cytokinesis 1 (PRC1) gene (NM_003981) at 15q26.1 (Figure 2c). This gene encodes the PRC1 protein, which is involved in cytokinesis and is a substrate for several cyclin-dependent kinases[22]. The PRC1 gene is down-regulated by the TP53 gene, and it is over-expressed in p53-defective cells[23]. Interestingly, the PRC1 gene is included in a five-gene expression signature that predicted prognosis among breast cancer patients in a recent study[24]. The expression level of the PRC1 gene was significantly higher in tumor tissue than in adjacent normal tissue (P = 4.62 × 10−30) among breast cancer cases included in TCGA (Supplementary Table 6). Our cis-eQTL analysis using TCGA data showed no association of rs2290203 with PRC1gene expression, but did reveal a correlation with expression of the RCCD1gene, which is 5,712 bp upstream of rs2290203. An eQTL analysis of human monocytes has also indicated that rs2290203 is a cis-eQTL for the RCCD1gene[16]. In our study, the rs2290203 risk allele (G) was associated with lower RCCD1 expression in both tumor (P = 3.6 × 10−4) and adjacent normal tissue (P = 0.007) (Supplementary Figure 2). However, these associations were no longer statistically significant after adjusting for the most significant cis-eQTL SNPs (rs4544218 for tumor tissue; rs59278520 for normal tissue), which are in strong LD with rs2290203 (Supplementary Figure 3). Variant rs4544218 was not associated with breast cancer risk (P = 0.8925), and rs59278520 was marginally associated with breast cancer risk (P = 0.0518) in the SBCGS-1, Stage 1 samples. The function of the RCCD1 gene is unknown.

SNP rs4951011 is located in intron 2 of the zinc finger CCCH domain-containing protein 11A (ZC3H11A) gene (NM_014827) at 1q32.1 (Figure 2a) and the 5′ untranslated region (UTR) of the zinc finger, BED-type containing 6 (ZBED6) gene (NM_001174108) (not shown in figure). The ZBED6 protein has recently been recognized as a novel transcription factor in placental mammals[25]. The function of the ZC3H11A gene is not clear. ChromHMM annotation using human mammary epithelial cell (HMEC) data from ENCODE suggests that rs4951011 may be located in a strong enhancer region marked by peaks of several active histone methylation modifications (H3K4me1, H3K4me3, H3K9ac, and H3K27ac). A search of RegulomeDB and HaploReg indicated that rs4951011 may be located in a predicted HNF1 motif and map to a DNase I hypersensitivity site in the MCF-7 cell line (Supplementary Table 7). Expression levels of the ZC3H11A gene were significantly higher in breast tumor tissue than in adjacent normal tissue (P = 0.0049) in TCGA data (Supplementary Table 6). Analyses using TCGA data revealed no evidence that this SNP or other SNPs correlated with it are cis-eQTLs for any genes in this locus. Recently, SNP rs4245739 in the MDM4 gene, ~752 kb downstream of rs4951011 (r = 0 in both ASN and CEU data), was associated with ER-negative breast cancer risk[5]. In our study, rs4245739 had a low MAF (0.03–0.05) and was not associated with breast cancer risk (P = 0.1861 in Stage 1).

In summary, our large GWAS conducted among East Asian women identified three new breast cancer susceptibility loci at 1q32.1, 5q14.3, and 15q26.1 and suggested a possible association with a fourth locus at 18q11.2. The associations of these loci with breast cancer risk may be mediated through cell-growth–control regulation, tumor-cell migration and invasion, or metastasis. Further studies of possible mechanisms through which these loci and genes are involved in breast tumorigenesis are warranted. Results from this study provide additional insights into the genetics and biology of breast cancer.

URLs

1000 Genomes Project, http://www.1000genomes.org/;

The Cancer Genome Atlas (TCGA), http://cancergenome.nih.gov/;

The DRIVE GAME-ON Consortium, http://epi.grants.cancer.gov/gameon/;

EIGENSTRAT, http://genepath.med.harvard.edu/~reich/EIGENSTRAT.htm;

ENCODE: http://www.genome.gov/10005107;

HaploReg, http://www.broadinstitute.org/mammals/haploreg/haploreg.php

HapMap Project, http://hapmap.ncbi.nlm.nih.gov/;

HumanExome Beadchip: http://genome.sph.umich.edu/wiki/Exome_Chip_Design;

LocusZoom, v1.1, http://csg.sph.umich.edu/locuszoom/;

MACH1.0, http://www.sph.umich.edu/csg/abecasis/MaCH/;

Mach2dat, http://genome.sph.umich.edu/wiki/Mach2dat:_Association_with_MACH_output;

METAL, http://www.sph.umich.edu/csg/abecasis/metal;

PLINK version 1.07, http://pngu.mgh.harvard.edu/~purcell/plink/;

R version 3.0.0, http://www.r-project.org/;

RegulomeDB, http://regulome.stanford.edu/;

SAS version 9.3, http://www.sas.com/;

UCSC Genome Browser, http://genome.ucsc.edu/;

ONLINE METHODS

Study populations

The ABCC comprises 22,780 cases and 24,181 controls from 14 studies (Supplementary Table 1), including 15,483 Chinese women, 18,367 Korean women, and 13,111 Japanese women. Chinese participants came from six studies based in: Shanghai [n = 12,219, the Shanghai Breast Cancer Study (SBCS), the Shanghai Breast Cancer Survival Study (SBCSS), the Shanghai Endometrial Cancer Study (SECS, controls only), and the Shanghai Women’s Health Study (SWHS)][6, 26–28], Taiwan (n = 2,131)[29], and Hong Kong (n = 1,133)[30]. Korean participants came from five studies: the Seoul Breast Cancer Study (SeBCS) (n = 6,179)[11], the Hwasun Cancer Epidemiology Study-Breast (HCES-Br) (n = 6,573)[31-33], the Korea Genome Epidemiology Study (KoGES; n = 3,209)[34], the Korean Hereditary Breast Cancer study (KOHBRA; n = 1,397)[35], and the Korean National Cancer Center (n = 1,009). Japanese participants came from three studies: the Biobank Japan Project (BBJ) (n = 11,021)[36], the Nagoya Study (n = 1,288)[37], and the Nagano Breast Cancer Study (n = 802)[38] (Supplementary Table 1). Detailed descriptions of these participating studies are presented in the Supplementary Note. The protocols for all participating studies were approved by their relevant Institutional Review Boards, and all participants of participating studies provided written informed consent. We estimated that our study had a statistical power of >80% to identify an association with an OR of 1.09 or above at P <5×10−8 for SNPs with a MAF as low as 0.25.

Genotyping Methods

Stage 1 genotyping

Stage 1 included two GWAS, in which 5,285 Chinese women and 4,777 Korean women were scanned primarily using the Affymetrix Genome-Wide Human SNP Array 6.0. Genotyping protocols for Stage I have been described elsewhere[6–9, 11]. In the Chinese GWAS (SBCGS-1), the initial 300 samples were genotyped using the Affymetrix GeneChip Mapping 500K Array Set. The remaining 4,985 samples were genotyped using the Affymetrix Genome-Wide Human SNP Array 6.0. After QC exclusions, the final dataset included 2,867 cases and 2,285 controls for 677,157 markers. For the Korean GWAS (SeBCS1), the Affymetrix Genome-Wide Human SNP Array 6.0 Array was also used[11]. After QC exclusions, the final dataset included 2,246 cases and 2,052 controls for 555,117 markers. Detailed descriptions of QC and criteria for sample and marker exclusion are presented in the Supplementary Note.

Stage 2 genotyping

Genotyping assays for 3,944 cases and 3,980 controls (SBCGS-2) in Stage 2 were completed using Illumina Infinium assays as the add-on content to the Illumina HumanExome Beadchip (see URLs). Genotype calling was carried out by using Illumina’s GenTrain version 2.0 clustering algorithm in GenomeStudio version 2011.1. Cluster boundaries were determined using study samples. Further QC procedures were conducted using PLINK (see URLs). Detailed descriptions of QC and criteria for sample and marker exclusion are presented in the Supplementary Note. Of the 4,598 SNPs selected, assays for 4,071 SNPs were successfully designed using Illumina Infinium assays. A total of 3,850 SNPs were successfully genotyped, and 3,678 SNPs were included in the analyses of 3,472 breast cancer cases and 3,595 controls.

Stage 3 genotyping

Genotyping assays for the 50 SNPs in Stage 3 were completed at the Vanderbilt Molecular Epidemiology Laboratory using the iPLEX Sequenom MassArray platform for 19,423 samples from the Taiwan, Hong Kong, HCES-Br, KOHBRA/KoGES, SeBCS2, Korea-NCC, Nagoya, and Nagano studies. Detailed descriptions of QC and criteria for sample and marker exclusion are presented in the Supplementary Note. For the BBJ1 study, the SNP data needed for the study were extracted from either genotyped (n = 8) or imputed (n = 14, mean RSQR = 0.96) data generated using the OmniExpress BeadChip. Breast cancer cases included in the BBJ2 study were genotyped using multiplex-PCR Invader assays. SNP data for the BBJ2 controls were extracted from data generated using the OmniExpress BeadChip.

We also selected 16 SNPs that showed a promising trend in the other studies included in Stage 3 for additional genotyping assays among 2,021 cases and 1,958 controls included in a case-control study conducted in Malaysia and Singapore that used the iPLEX Sequenom MassArray platform at the Cancer Research Initiatives Foundation, Sime Darby Medical Centre, Malaysia. However, because of a potential concern about genetic admixture revealed in our previous study[3] and an unusual pattern of associations observed in these studies (Supplementary Table 8), we did not include these samples in the final analysis.

Statistical analyses

PLINK version 1.07 (see URLs)[39] was used to analyze the genome-wide data obtained in Stage 1. To evaluate the population structure in the Chinese GWAS (SBCGS-1), we performed principal component analyses using the EIGENSTRAT software[40] in a set of approximately 6,000 independent SNPs that met the following criteria: (1) a neighboring distance > 200 kb from the SNP of interest, (2) a MAF > 0.2, (3) r2 < 0.1, and (4) a genotype call rate > 99%. The inflation factor (λ) was estimated to be 1.0426. Similar analyses were performed for the Korean GWAS (SeBCS1) and yielded a λ of 1.0431[11]. We also rescaled the inflation statistic to an equivalent value for a study including 1,000 cases and 1,000 controls (λ1000) using the formula: λ1000 = 1 + 500 × (λ − 1) × (1/Ncases + 1/Ncontrols)[41]. The λ1000 was 1.02 for both studies included in Stage 1 and 1.01 in the meta-analysis of both studies. These data suggest that any population substructure, if present, should not have any appreciable effect on the results. ORs associated with each SNP and 95% confidence intervals (CIs) were estimated using logistic regression implemented in PLINK with adjustment for age and the first two principal components.

We used the program MACH 1.0 (see URLs)[12] to impute genotypes for autosomal SNPs (n = 2,416,663) that were present in CHB (Han Chinese in Beijing, China) and JPT (Japanese in Tokyo, Japan) HapMap Phase II release 22 data for samples included the Chinese and Korean GWAS. Only SNPs with high imputation quality score (RSQR > 0.50) and a MAF > 0.05 in these two GWAS were included in the analyses. Dosage data for imputed SNPs for samples in each GWAS were analyzed using the program Mach2dat (see URLs)[12]. Associations between genotype dosage (0, 1, and 2) of the effect allele and breast cancer risk were assessed using logistic regression models after adjusting for age and the first two principal components. ORs associated with each SNP and 95% CIs were estimated under a log-additive model. We also used SAS version 9.3 (see URLs) to analyze genotype data, which yielded results virtually identical to those generated from dosage data using Mach2dat. We obtained summary ORs and 95% CIs of SNPs from the two GWAS by using METAL software (see URLs)[13] to run a fixed-effects inverse-variance meta-analysis.

Individual data were obtained from all studies except for the two BioBank Japan studies (BBJ1 and BBJ2). Case-control differences in selected demographic characteristics and major risk factors were evaluated using t-tests (for continuous variables) and Chi-square tests (for categorical variables). Summary associations between SNPs and breast cancer risk were generated on the basis of a fixed-effects inverse-variance meta-analysis conducted using METAL software[13]. Analyses stratified by ethnicity and estrogen receptor (ER) status were also carried out. Heterogeneity across studies, among ethnicities, and according to ER status was assessed with a Cochran’s Q test. P values of < 5 × 10−8 in the combined analysis were considered statistically significant.

We assessed associations of breast cancer risk with the three newly identified risk variants among European-ancestry women in collaboration with the DRIVE GAME-ON Consortium (see URLs). Included in this analysis were data from 16,003 cases and 41,335 controls recruited in twelve studies. Genome-wide scan data from these studies were imputed and meta-analyzed, and summary data are presented herein.

We generated forest plots by using R version 3.0.0 (see URLs). Regional association plots were drawn using the website-based tool LocusZoom, version 1.1 (see URLs)[42]. LD matrices used in this study were reported based on HapMap release 22 data. All genomic references are based on NCBI Build 36.

Functional annotation

Functional annotation was performed using data from the ENCODE project (see URLs) accessed through the UCSC Genome Browser (see URLs). Enhancer and transcription elongation regions were predicted in HMEC by using ChromHMM annotation. DNase I hypersensitive areas, transcription factor binding sites, and miRNA binding sites were evaluated in all cell types, including breast cancer cells, included in ENCODE. RegulomeDB (see URLs)[17], a database that annotates SNPs with known and predicted regulatory elements in the intergenic regions of the human genome using data from GEO, the ENCODE project, and published literature, was also used to predict the possibility of transcription factor binding sites and DNase I hypersensitivity. In addition, we performed functional annotation using HaploReg v2 (see URLs)[18], a tool for exploring annotations of the noncoding genome at variants on haplotype blocks.

eQTL analysis

We used the TCGA breast cancer data (Supplementary Note) to perform an eQTL analysis for normal and tumor tissue samples separately. Detailed descriptions of eQTL analysis are presented in the Supplementary Note. We focused only on the SNPs and genes located within the 1Mb regions flanking the three newly identified risk loci to identify cis-eQTLs. A significance threshold P value of < 0.01 was used to determine candidate cis-eQTLs.

Differential gene expression analysis

To identify differentially expressed genes located in the three identified risk loci, we analyzed data from a total of 87 pairs of breast tumor-normal tissue samples included in TCGA (Supplementary Note). Detailed descriptions of differential gene expression analysis are presented in the Supplementary Note.