Trans-ethnic predicted expression genome-wide association analysis identifies a gene for estrogen receptor-negative breast cancer.

Genome-wide association studies (GWAS) have identified more than 90 susceptibility loci for breast cancer, but the underlying biology of those associations needs to be further elucidated. More genetic factors for breast cancer are yet to be identified but sample size constraints preclude the identification of individual genetic variants with weak effects using traditional GWAS methods. To address this challenge, we utilized a gene-level expression-based method, implemented in the MetaXcan software, to predict gene expression levels for 11,536 genes using expression quantitative trait loci and examine the genetically-predicted expression of specific genes for association with overall breast cancer risk and estrogen receptor (ER)-negative breast cancer risk. Using GWAS datasets from a Challenge launched by National Cancer Institute, we identified TP53INP2 (tumor protein p53-inducible nuclear protein 2) at 20q11.22 to be significantly associated with ER-negative breast cancer (Z = -5.013, p = 5.35×10-7, Bonferroni threshold = 4.33×10-6). The association was consistent across four GWAS datasets, representing European, African and Asian ancestry populations. There are 6 single nucleotide polymorphisms (SNPs) included in the prediction of TP53INP2 expression and five of them were associated with estrogen-receptor negative breast cancer, although none of the SNP-level associations reached genome-wide significance. We conducted a replication study using a dataset outside of the Challenge, and found the association between TP53INP2 and ER-negative breast cancer was significant (p = 5.07x10-3). Expression of HP (16q22.2) showed a suggestive association with ER-negative breast cancer in the discovery phase (Z = 4.30, p = 1.70x10-5) although the association was not significant after Bonferroni adjustment. Of the 249 genes that are 250 kb within known breast cancer susceptibility loci identified from previous GWAS, 20 genes (8.0%) were statistically significant associated with ER-negative breast cancer (p<0.05), compared to 582 (5.2%) of 11,287 genes that are not close to previous GWAS loci. This study demonstrated that expression-based gene mapping is a promising approach for identifying cancer susceptibility genes.

Introduction

Breast cancer is the most common cancer in women in the United States and in the world [1]. It is a heterogeneous disease and the two main subgroups of breast cancer are estrogen receptor (ER)-positive and ER-negative cancer. Genome-wide association studies (GWAS) have identified more than 90 susceptibility loci for breast cancer [2-20], with only a few loci specific for ER-negative breast cancer [3,15,17]. Susceptibility loci for ER-positive loci are often the same as loci for overall breast cancer risk because most of breast cancers are ER-positive, especially in women of European or Asian ancestry [2,4,19].

Women of African ancestry are more likely to be diagnosed with ER-negative breast cancer compared to women of non-African ancestry [21-23]. To date, breast cancer GWAS have been conducted primarily in populations of European ancestry. The difference in linkage disequilibrium (LD) patterns and allele frequencies across ancestry groups may explain the apparent inconsistencies in GWAS findings from studies of women of European ancestry as compared to studies of women of African ancestry [24,25]. The strength and the direction of the association between causal variants and disease are expected to be consistent across populations, and thus cross-population validation provides further evidence of causation. In addition, trans-ancestry analysis could identify novel breast cancer susceptibility variants [26].

The variants discovered by previous GWAS along with previously known high-penetrance genes explain only a modest proportion of the heritability of breast cancer [2]. More genetic factors for breast cancer are yet to be identified, but power for discovery of new loci is limited by the sample size of existing GWASs. Moreover, the biologic significance of the variants identified by GWAS and the genes on which they act, are often unknown. Single nucleoid polymorphisms (SNPs) associated with disease traits are more likely to be expression quantitative trait loci (eQTLs) [27], and regulatory variants can explain a large proportion of disease heritability [28]. Therefore, genes regulated by eQTLs can be used as an enrichment analysis unit to identify more genetic risk factors for breast cancer. Recently, gene-based approaches using eQTL information, such as PrediXcan, have been proposed, which can reduce the multiple testing burden in genome-wide analyses and have been used to identify novel genes for autoimmune diseases [29]. PrediXcan uses individual-level data to estimate the correlation between genetically predicted levels of gene expression and human traits to prioritize causal genes. MetaXcan computes the same correlation as PrediXcan, but does so using summary statistics from GWAS, which are much more readily accessible than individual level data [30].

To identify novel genes involved in breast cancer susceptibility, we utilized a gene-level expression-based association method, implemented in the MetaXcan software [30], to infer gene expression levels using summary statistics from five GWASs. We used an additive prediction model of gene-expression levels trained in Depression Genes and Network (DGN) data [31] and examined the predicted expression of specific genes for association with overall breast cancer risk and estrogen receptor-negative breast cancer risk. The GWAS datasets were made available in dbGaP (https://www.ncbi.nlm.nih.gov/gap) through “Up For A Challenge (U4C)–Stimulating Innovation in Breast Cancer Genetic Epidemiology” launched by the National Cancer Institute. The DGN data included RNA sequencing data from whole blood of 922 genotyped individuals (463 cases of major depressive disorder and 459 controls), all of European ancestry. These individuals consisted of 274 males and 648 females with ages ranged from 21 to 60.

Results

Using logistic regression, we first conducted SNP-level GWAS analysis for overall breast cancer risk among 8605 breast cancer cases and 8095 controls, and for ER-negative breast cancer risk among 3879 cases and 10213 controls. The analyses were performed for each of the five GWAS datasets separately and summary statistics including log odds ratios and standard errors were generated. These summary statistics for each dataset were input to the software MetaXcan [30] to perform genome-wide gene-level expression association tests for 11,536 genes. Then, we performed meta-analysis of the results from individual MetaXcan analyses. Quantile-quantile plots of P-values from the meta-analysis showed little inflation (). For overall breast cancer risk, there was no gene with a P-value that deviated from the null distribution (), but for ER-negative breast cancer risk analysis, there were several genes with P-values smaller than expected, including TP53INP2, HP, and DHODH ().

Quantile-quantile plot of gene-based association P values for (a) overall and (b) estrogen receptor negative breast cancer. Red line shows the null distribution of P values.

lists the top genes with P-values less than 10−3 in the analyses of association between predicted gene expressions and overall breast cancer risk. The sign of Z score indicates the direction of association between genetically-predicted expression and breast cancer risk. None of the genes reached genome-wide significance when a Bonferroni threshold (α = 4.33x10-6) was used.

Of the 249 genes that are 250 kb within known susceptibility loci identified from previous breast cancer GWAS [2-4,17,32], 12 genes (4.8%) were statistically significant associated with overall breast cancer risk at nominal significance level of 0.05, compared to 497 (4.4%) of 11,287 genes that are not close to previous GWAS loci (P for enrichment = 0.75).

lists the genes with P-values less than 10−3 in the ER-negative breast cancer analysis. TP53INP2 was the top gene (P = 5.35x10-7), which surpassed the Bonferroni-corrected p-value threshold (α = 4.33x10-6). The false discovery rate for TP53INP2 was 0.0062. Higher genetically-predicted TP53INP2 expression was associated with lower risk of ER-negative breast cancer. The gene with the second smallest P-value was HP, which had p-value of 1.70x10-5, close to but not significant after Bonferroni correction. The false discovery rate for the HP gene was 0.098. For the HP gene, higher expression was associated with higher risk of ER-negative breast cancer. Both genes are novel and no previous studies have found association between these two genes and breast cancer risk.

Of the 249 genes that are 250 kb within known breast cancer susceptibility loci identified from previous GWAS, 20 genes (8.0%) were statistically significant associated with ER-negative breast cancer (p<0.05), compared to 582 (5.2%) of 11,287 genes that are not close to previous GWAS loci (P for enrichment = 0.044), suggesting a moderate enrichment for genes close to known susceptibility loci.

There were six SNPs included in the prediction of the expression of the TP53INP2 gene, from 367 kb upstream to 159 kb downstream of the gene (). Five of the six SNPs (except for rs8116198) were associated with overall breast cancer risk and ER-negative breast cancer risk (at the nominal level of α = 0.05), and the effects were consistently across studies (none of the heterogeneity tests were significant). These associations were more significant for ER-negative breast cancer risk (p values ranging from 5.0x10-4 to 1.8x10-6) than for overall breast cancer risk (7.0x10-4 to 1.4x10-4). None of the SNP-level associations reached traditional genome-wide significance, thus they have not been reported in previous GWAS publications. However, our study showed the aggregate effects of these SNPs were significantly associated with ER-negative breast cancer after Bonferroni correction. We noticed that one of the six SNPs, rs8116198, is monomorphic in the SBCGS data. Therefore, when MetaXcan was applied to the SBCGS data, the prediction of TP53INP2 expression was based on only five SNPs. To make our results more robust to missing and low quality genotypes, in the DGN prediction model, we used elastic net with 0.5 as the mixing parameter, which sets the degree of mixing between ridge regression and LASSO. In addition, the SNPs in the prediction were not necessarily causal but could be in LD with the causal SNPs.

shows positions of the 6 eQTL SNPs for TP53INP2 in the cytoband 20q11.22. Interestingly, there are several other genes in this region that were associated with ER-negative breast cancer, including MAP1LC3A, ITCH, and TRPC4AP ( and ). The 6 SNPs are located either in enhancer elements or in promotor regions (). The promotor/enhancer features of 4 SNPs were found in human mammary epithelial cells (HMEC) and breast variant human mammary epithelial cells (HMEC.35), and the enrichment was statistically significant for both cell types (both p<0.03).

The 20q11.22 locus spanned for expression quantitative trait loci of TN53INT2, and analysis of regulation enhancer with data from ENCODE through UCSC Genome Browser, including transcription factor binding sites and human mammary epithelial cells (HMEC) histone modification marks.

Chromosomal coordinates are in NCBI build 37.

There were 20 SNPs included in the prediction of the expression of the HP gene (). Thirteen of the 20 SNPs were associated with overall breast cancer risk and 17 were associated with the risk of ER-negative breast cancer (at the nominal level of α = 0.05), quite consistently across populations (none of the heterogeneity tests were significant). The strengths of their associations were stronger for ER-negative breast cancer risk than for overall breast cancer risk. Interestingly, none of the associations for individual SNPs reached genome-wide significance, thus they have not been reported in previous GWAS publications.

We used summary results from GAME-ON GWAS (http://gameon.dfci.harvard.edu) to replicate our study findings from the U4C. All the six eQTLs for the TP53INP2 gene were available in GAME-ON (). Five of the six SNPs that were associated with ER-negative breast cancer in the discovery phase (using U4C datasets) were all statistically significant in GAME-ON at the nominal 0.05 significance level. Gene-level test of TP53INP2 from MetaXcan gave a Z-score of -2.803 (p = 5.1×10−3) for ER-negative breast cancer in GAME-ON. The gene-level test for overall breast cancer risk was not significant in GAME-ON (Z-score = -1.627, p = 0.10). Because the GAME-ON ER-negative data included the BPC3 dataset, in order to show the independent replication, we tested association in the U4C ER-negative data excluding BPC3, and found the Z-score for the TP53INP2 gene was -4.127 (p = 3.67×10−5).

For the HP gene, the direction of association for 19 SNPs (out of 20) were consistent between U4C and GAME-ON for ER-negative breast cancer risk, but only 2 SNPs were statistically significant at nominal 0.05 level in GAME-ON (). None of the SNPs were significantly associated with overall breast cancer risk in GAME-ON. In the gene-based analysis using GAME-ON data, the Z-score for overall breast cancer risk was 1.769 (p = 0.077) and the Z-score for ER-negative breast cancer risk was 2.02 (p = 0.043). In addition, we tested this association in the U4C ER-negative data excluding BPC3, and found the Z-score for the HP gene was 2.81 (p = 5.1×10−3).

Discussion

In this gene-level expression-based genome-wide association analysis of five breast cancer GWAS datasets composed of individuals of diverse ancestry, we identified TP53INP2 (20q11.22) as gene with genetically-determined expression that is associated with ER-negative breast cancer. The gene-based analysis of aggregated eQTLs for a particular gene as an analysis unit can reduce the burden of multiple testing and provide a direction of association between expression of a specific gene and disease risk. We found that increased expression of TP53INP2 expression in whole blood was associated with a decrease in ER-negative breast cancer risk. In addition, we identified the HP gene in the 16q22.2 regions to have expression levels that are positively associated with ER- negative breast cancer.

The TP53INP2 gene (tumor protein p53-inducible nuclear protein 2) is 9150 base pairs long and codes for a 220 amino acid protein, which is a dual regulator of transcription and autophagy and is required for autophagosome formation and processing. One experimental study showed that overexpression of TP53INP2 severely attenuated proliferative and invasive capacity of melanoma cells, via p53 signaling and lysosomal pathways [34]. This inverse correlation between TP53INP2 expression and cancer proliferation is consistent with our finding that TP53INP2 expression inversely correlated with breast cancer risk. P53 is a transcription factor for TP53INP2, and TP53 plays an important role in development of multiple cancers. Germline TP53 mutations cause Li-Fraumeni syndrome, characterized as a cluster of cancers including breast cancer [35]. Somatic TP53 mutation is a common event in ER-negative breast cancer [36]. As a downstream gene of p53, TP53INP2 may affect breast cancer risk through p53 signaling pathway. Also, known as DOR (diabetes- and obesity-regulated gene), TP53INP2 has been linked to obesity and diabetes [37]. TP53INP2 is also associated with triglycerides and cholesterol level. One experimental study found that dietary fat content influenced the expression of TP53INP2 expression in adipose and muscle tissues of mice [38]. This gene has been proposed to serve as a diagnostic biomarker for papillary thyroid carcinoma [39] but no study has linked its expression to cancer risk. Obesity has been convincingly correlated with breast cancer risk in numerous studies, although the relationship is complex and involves additional modifying factors [40,41]. Obesity has been associated with excess risk for breast cancer among postmenopausal women [42-46], while in pre-menopausal women, obesity was associated with decreased breast cancer risk [40,43,47-49]. However, the underlying mechanisms for this association are still not fully understood. The identification of TP53INP2/DOR as breast cancer-related gene could provide novel insight on the mechanism for obesity-breast cancer relationship.

In the 20q11.22 region, several other genes including MAP1LC3A, ITCH, and TRPC4AP were associated with ER-negative breast cancer risk. MAP1LC3A codes for a protein that is important in the autophagy process, and was found to be expressed at higher level in breast cancer tissues than in normal tissues [50]. E3 ubiquitin ligase ITCH plays a role in erythroid and lymphoid cell differentiation and immune response regulation, and ITCH was found to be important in the cross-talk between the Wnt and Hippo pathways in breast cancer development [51]. TRPC4AP is involved in Ca2+ signaling and is part of the ubiquitin ligase complex [52,53]. It is unclear which of these genes (or their interactions) play a role in breast cancer development, but the 20q11.22 locus is worthy of further investigation. Three of the six SNPs for TP53INP2 (rs6060047, rs11546155, and rs1205339) are also shared by the genes MAP1LC3A and TRPC4AP. It is possible that the associations in these three genes are partly generated by the overlapped SNPs, which contribute to predicted expression levels of the three genes and, possibly, to the enrichment observed at this locus.

The HP gene (16q22.2) is 6,491 base pairs long and codes for a 406 amino acid preprotein, which codes haptoglobin. Haptoglobin binds to hemoglobin to prevent iron loss during hemolysis. There are two allelic forms, Hp1 (83 residues) and Hp2 (142 residues), which determine 3 major phenotypes [54]. Haptoglobin genotype has been linked to cardiocerebral outcomes among diabetic patients [55]. A small study found haptoglobin phenotypic polymorphism was associated with familial breast cancer [56], but no studies have reported on the relationship between SNPs in this gene and breast cancer risk. Further larger studies could investigate the relationship between major HP genotype/phenotype (HP1-1, HP1-2, and HP2-2) and breast cancer risk.

The present study has several strengths, including its large sample size, diverse ancestry groups, a cross-replication approach, and a novel gene expression-based analysis method. The gene-level analysis method can combine eQTL SNPs in a biologically informative way to assess relationships between predicated gene expression and disease risk. Compared to SNP-based analysis, the gene-based analysis can gain power by reducing the multiple testing burden by about 100-fold and using external information on correlation between gene expression and SNPs from reference samples. In addition, this approach enables the detection of individual SNPs with weak effects on disease risk by leveraging combined effects of multiple SNPs on gene expression. For example, none of SNPs for TP53INP2 reached traditional genome-wide significance, but their aggregated effect via TP53INP2 expression was genome-wide significant. The gene-based method (MetaXcan) that we employed is an extension of the gene expression-based method (PrediXcan) [29] and allows the use of SNP-level summary statistics without the need to access individual-level genotype data [30]. The MetaXcan method has been shown to produce PrediXcan results accurately, and it is robust to ancestry mismatches between studies and reference/training populations [30]. With this property, we were able to use summary statistics from the GAME-ON consortium for external replication.

Several limitations should be considered when interpreting the study findings. The gene expression-based association method relies on accurate prediction of gene transcript level from genotypes, i.e. identification of eQTLs, but eQTL identification depends on sample size of eQTL studies as well as tissue types. In the current study, we used the transcriptome prediction model that was developed using 922 RNA-seq samples from whole blood and genotype data [31]. Although it has been shown that models developed in whole blood were still useful for understanding diseases that affect other primary tissues [29], we expect there to be a loss of power when studying non-blood diseases using whole blood eQTL data. As a sensitivity analysis, we performed the MetaXcan analysis using the prediction model from breast tissues of 183 donors of multiple ethnicities (http://www.gtexportal.org). Only 4,308 genes had breast tissue specific eQTLs, and no eQTL was available for TP53INP2, perhaps due to the small sample size. We found that DHODH (P = 3.61×10−5), ITCH (P = 1.23×10−4), and TRPC4AP (P = 7.7x10-4) were among the top genes associated with ER-negative breast cancer risk, and TRPC4AP (P = 1.68×10−5) and DHODH (P = 1.12×10−4) among the top genes associated with the overall breast cancer risk using breast tissue eQTLs. In the enrichment analysis, we found that 7 (8.2%) out of 85 genes that are close to known breast cancer susceptibility loci identified in previous GWAS were associated with ER-negative breast cancer and 6 (7.1%) genes were associated with overall breast cancer risk; by contrast, of the 4223 genes away from previous GWAS loci, 199 (4.7%) genes were associated with ER-negative breast cancer and 212 (5.0%) genes were associated with overall breast cancer risk. Here, we have to consider the balance between tissue relevance and sample size in eQTL studies. Further investigations based on large, reliable eQTL datasets are desirable. In future studies, we will seek out larger samples of multi-ethnic breast tissue as training data to construct improved prediction models of gene expression and further investigate trans-ethnic associations for breast cancer.

In conclusion, our study identified TP53INP2 and several other genes in the 20q11.22 region as potential susceptibility genes for ER-negative breast cancer using a novel gene-based analysis method that incorporates genetically determined gene expression. We demonstrated this gene-based method increases statistical power and may be helpful in searching for causal variants. Future studies need to determine whether the TP53INP2 gene is a true susceptibility gene for breast cancer and what are the underlying mechanisms for its association with ER-negative breast cancer.

Materials and methods

Study samples

The study was approved by the Institutional Review Board of the University of Chicago. The Epidemiology and Genomic Research Program within the National Cancer Institute launched a Challenge at the end of 2015 to inspire novel cross-disciplinary approaches to more fully decipher the genomic basis of breast cancer, called "Up For A Challenge (U4C)–Stimulating Innovation in Breast Cancer Genetic Epidemiology”. Several data sets were gathered and made available for use in dbGap (https://www.ncbi.nlm.nih.gov/gap). Our study has two phases; the discovery phase included five U4C GWAS datasets (). Here, we refer them collectively as “U4C” data. These data were collected from three distinct ancestry groups. The BPC3 [16,18] and CGEMS study [15,20] were conducted in women of European ancestry. The ROOT [17] and AABC study [57] consisted of women of African ancestry. The SBCGS study was conducted in Chinese population [19]. For the analysis of overall breast cancer risk, we used four GWAS datasets: AABC, CGEMS, ROOT, and SBCGS. For the analysis of ER-negative breast cancer risk, we used datasets from AABC, BPC3, ROOT, and SBCGS. All these dbGap datasets included imputed genotype data that were inferred based on reference haplotypes from the 1000 Genomes project.

In the replication phase, we used the summary results from the meta-analysis of 11 breast cancer GWASs in the GAME-ON consortium (http://gameon.dfci.harvard.edu). All participants were of European ancestry. The overall breast cancer analysis included 16,003 cases and 41,335 controls from 11 GWAS studies; The ER-negative breast cancer analysis included 4939 cases and 13128 controls from 7 GWAS studies. The dataset from one study (BPC3; all ER-negative cases) in GAME-ON consortium was also included the U4C datasets. Because only meta-analysis results were available from GAME-ON, we removed the BPC3 data from “U4C” dataset when we compared replication performance to avoid duplicate counting.

Statistical analysis

Our gene level expression-based association analysis consists of three main steps. First, we conducted SNP-level genome-wide association tests and calculated summary statistics such log odds ratios and their standard errors. We used logistic regression model adjusting for eigenvectors from the principal component analysis and related covariates such as age. Genotypes were coded by an additive genetic model. Eigenvectors in principal component analysis were calculated using the method smartPCA, which is implemented in the software EIGENSOFT version 6.0.1 [58]. For the ROOT dataset, we adjusted for age, study sites, and the top 4 eigenvectors. For the AABC dataset, we adjusted for age, study sites, and top 10 eigenvectors. For CGEMS and SBCGS, we adjusted for age and the top three or two eigenvectors, respectively. The number of eigenvectors we adjusted for was chosen according to published papers from these GWASs [17,57], as well as their association with case-control status. The logistic regression models were fit using software Mach2dat (http://www.unc.edu/~yunmli/software.html) or SNPtest [59], depending on format of the datasets; the Mach2dat software was used for CGEMS and SBCGS and SNPtest was used for ROOT and AABC. For BPC3, the GWAS summary statistics for ER-negative breast cancer have been pre-calculated in the dbGap release, so we used them directly.

Second, we applied the gene level association method, MetaXcan [30] (https://github.com/hakyimlab/MetaXcan), to each of the datasets listed in . MetaXcan is an extension of the method PrediXcan [29], which uses an additive genetic model to estimate the component of gene expression determined by an individual’s genetic profile and then identifies likely causal genes by computing the correlations between genetically predicted gene expression levels and disease phenotypes. MetaXcan infers the results of PrediXcan using summary statistics from GWAS, which are much more readily accessible than individual level data. In our study, as input files for MetaXcan, we used summary statistics from SNP-based analysis of each dataset obtained in step one. In addition, we used the whole blood genetic prediction model of transcriptome levels trained in the DGN data [31], which can be downloaded from http://predictdb.hakyimlab.orghttps://s3.amazonaws.com/predictdb/DGN-HapMap-2015/. The DGN data provides a large reference sample of 922 individuals with both genome-wide genotype data and RNA sequencing data. The model trained in the DGN data can be useful in estimating gene expression levels and has been successfully applied to the Wellcome Trust Case Control Consortium (WTCCC) data in identifying genes associated with five complex diseases [29]. The DGN prediction model includes a) weights for predicting gene expression using genotypes and b) covariance of the SNPs that takes into account linkage disequilibrium. We tested the association between predicted expression levels of 11,536 genes for each of the two phenotypes, overall and ER-negative breast cancer risk, using the MetaXcan software. To construct the prediction model of expression levels using the DGN data, MetaXcan used SNPs with minor allele frequencies (MAFs) >0.05. When MetaXcan was applied to the breast cancer GWAS data, only SNPs with MAFs >0.05 were used. We also looked up genes within 250 kb of the 93 breast cancer susceptibility loci identified in previous GWAS [2-4,17,32].

Third, we conducted meta-analysis to combine results from MetaXcan analyses for different datasets. The method described by Willer et al. with sample size as meta-analysis weight [60] was used. We also conducted SNP-level meta-analysis using a fixed effect model, as implemented in the software METAL (http://genome.sph.umich.edu/wiki/METAL). False discovery rates were calculated using the Benjamini and Hochberg method [61].

For genes identified in the discovery phase using the U4C datasets, we conducted replication analysis using GAME-ON summary results using the same methods described above. For each top variant and gene identified in this study, we used HaploReg [33] and USCS Genome Browser to explore functional annotations of noncoding variants. Chromatin states (promoters and enhancers), variant effect on regulatory motifs, and protein binding sites were assessed from available data from the ENCODE [62] and Roadmap Epigenomics Consortium [63]. Data from normal mammary epithelial cells (HMEC, MYO, vMHEC) were emphasized.

HP-related SNPs and their association with breast cancer risk.

(DOCX)

Click here for additional data file.

GAME-ON replication for SNPs related to the HP gene.

(DOCX)

Click here for additional data file.

Table 1

Top genes with P-values < 10−3 in analyses of association between predicted gene expressions and overall breast cancer risk*.

Gene	Cytoband	SNPs in predictor	AABC		CGEMS		ROOT		SBCGS		Total
			Z score	P value	Z score	P value	Z score	P value	Z score	P value	Z score	P value	FDR
TP53INP2	20q11.22	6	-2.536	1.12E-02	-0.953	3.41E-01	-3.023	2.50E-03	-1.683	9.23E-02	-4.180	2.91E-05	0.34
BAG3	10q26.11	18	-2.109	3.49E-02	-1.145	2.52E-01	-3.003	2.67E-03	-1.074	2.83E-01	-3.660	2.52E-04	0.77
POLN	4p16.3	39	-2.291	2.20E-02	-2.614	8.96E-03	-0.942	3.46E-01	-1.624	1.04E-01	-3.644	2.68E-04	0.77
WDR37	10p15.3	9	-1.637	1.02E-01	-0.747	4.55E-01	-0.292	7.70E-01	-4.144	3.42E-05	-3.629	2.84E-04	0.77
TTLL5	14q24.3	26	2.717	6.59E-03	2.087	3.69E-02	1.189	2.35E-01	1.206	2.28E-01	3.588	3.33E-04	0.77
HP	16q22.2	19	2.424	1.53E-02	1.961	4.99E-02	1.598	1.10E-01	1.147	2.52E-01	3.529	4.18E-04	0.77
VTI1B	14q24.1	1	1.433	1.52E-01	2.790	5.28E-03	0.902	3.67E-01	2.151	3.15E-02	3.471	5.18E-04	0.77
HLA-DMA	6p21.32	30	-2.001	4.54E-02	-0.756	4.50E-01	-1.603	1.09E-01	-2.293	2.19E-02	-3.456	5.47E-04	0.77
MYOM2	8p23.3	109	2.338	1.94E-02	-0.051	9.59E-01	2.153	3.14E-02	1.955	5.06E-02	3.430	6.04E-04	0.77
MYO9B	19p13.11	6	1.643	1.00E-01	0.887	3.75E-01	1.473	1.41E-01	2.549	1.08E-02	3.373	7.44E-04	0.81
ZNF202	11q24.1	15	2.214	2.69E-02	1.675	9.39E-02	0.003	9.98E-01	2.644	8.20E-03	3.363	7.71E-04	0.81

*Bonferroni threshold = 4.33×10−6.

FDR, false discovery rate.

Table 2

Top genes with P-values < 10−3 in analyses of association between predicted gene expressions and ER-negative breast cancer risk*.

Gene	Cytoband	SNPs in predictor	AABC		BPC3		ROOT		SBCGS		Total
			Z score	P value	Z score	P value	Z score	P value	Z score	P value	Z score	P value	FDR
TP53INP2	20q11.22	6	-3.708	2.09E-04	-2.919	3.51E-03	-2.703	6.87E-03	-0.417	6.77E-01	-5.013	5.35E-07	0.0062
HP	16q22.2	20	1.424	1.54E-01	3.302	9.61E-04	1.851	6.41E-02	1.749	8.03E-02	4.300	1.70E-05	0.098
DHODH	16q22.2	58	-1.121	2.62E-01	-4.700	2.61E-06	1.020	3.08E-01	-1.859	6.31E-02	-4.119	3.80E-05	0.15
YJEFN3	19p13.11	20	-2.650	8.05E-03	-2.797	5.16E-03	0.154	8.78E-01	-1.549	1.21E-01	-3.810	1.39E-04	0.34
MAP1LC3A	20q11.22	49	-2.077	3.78E-02	-2.922	3.48E-03	-1.751	7.99E-02	-0.157	8.75E-01	-3.734	1.88E-04	0.34
DPY19L1	7p14.2	24	2.188	2.87E-02	3.035	2.41E-03	0.672	5.01E-01	0.791	4.29E-01	3.731	1.91E-04	0.34
GCOM1	15q21.3	75	-1.841	6.56E-02	-3.295	9.85E-04	-0.854	3.93E-01	-0.525	5.99E-01	-3.689	2.25E-04	0.34
AMOTL1	11q21	14	2.155	3.12E-02	2.118	3.42E-02	0.448	6.54E-01	2.509	1.21E-02	3.675	2.38E-04	0.34
ITCH	20q11.22	12	-1.597	1.10E-01	-3.861	1.13E-04	0.203	8.39E-01	-0.318	7.51E-01	-3.494	4.77E-04	0.61
TRPC4AP	20q11.22	26	2.385	1.71E-02	1.899	5.76E-02	2.536	1.12E-02	0.127	8.99E-01	3.466	5.28E-04	0.61
SNX24	5q23.2	3	-0.902	3.67E-01	-2.235	2.54E-02	-1.612	1.07E-01	-2.022	4.32E-02	-3.327	8.77E-04	0.91

*Bonferroni threshold = 4.33×10−6

FDR, false discovery rate

Table 3

TP53INP2-related SNPs and their association with breast cancer risk.

	Pos. at chr20/from TP53INP2*			Overall		ER-negative
SNP		Test/ref allele	Study	OR (95% CI)	P	OR (95% CI)	P
rs1205339	32,924,967	G/A	BPC3			1.17 (1.05–1.31)	5.9E-03
	-367,127		CGEMS	1.05 (0.91–1.23)	0.49
			AABC	1.14 (1.03–1.27)	0.013	1.29 (1.12–1.50)	5.6E-04
			ROOT	1.19 (1.04–1.36)	0.012	1.39 (1.10–1.76)	5.6E-03
			SBCGS	1.07 (0.97–1.18)	0.18	1.06 (0.89–1.26)	0.51
			meta	1.11 (1.05–1.18)	4.2E-04	1.20 (1.11–1.29)	1.8E-06
rs4911154	32,996,101	A/G	BPC3			1.16 (1.04–1.30)	0.01
	-295,993		CGEMS	1.04 (0.89–1.21)	0.65
			AABC	1.15 (1.03–1.28)	0.014	1.32 (1.14–1.54)	2.9E-04
			ROOT	1.23 (1.07–1.42)	3.5E-03	1.39 (1.09–1.78)	8.3E-03
			SBCGS	1.09 (0.98–1.21)	0.11	1.06 (0.88–1.28)	0.54
			meta	1.13 (1.06–1.20)	1.6E-04	1.20 (1.11–1.30)	2.6E-06
rs8116198†	33,114,201	G/A	BPC3			0.92 (0.80–1.05)	0.21
	-177,893		CGEMS	0.94 (0.78–1.13)	0.5
			AABC	1.04 (0.82–1.33)	0.73	1.11 (0.79–1.57)	0.54
			ROOT	0.67 (0.45–1.00)	0.052	0.87 (0.49–1.54)	0.63
			meta	0.93 (0.81–1.07)	0.33	0.92 (0.80–1.04)	0.30
rs6058107	33,288,546	C/T	BPC3			0.87 (0.78–0.97)	8.7E-03
	-3,548		CGEMS	0.92 (0.80–1.06)	0.27
			AABC	0.91 (0.83–1.01)	0.072	0.84 (0.73–0.96)	0.014
			ROOT	0.83 (0.73–0.94)	4.00E-03	0.80 (0.64–0.99)	0.041
			SBCGS	0.92 (0.84–1.00)	0.057	1.03 (0.88–1.20)	0.71
			meta	0.90 (0.85–0.95)	1.4E-04	0.90 (0.83–0.97)	5.0E-04
rs6060047	33,367,400	G/T	BPC3			0.87 (0.77–0.97)	0.017
	75,306		CGEMS	0.91 (0.78–1.07)	0.27		.
			AABC	0.88 (0.79–0.98)	0.016	0.75 (0.65–0.86)	7.5E-05
			ROOT	0.83 (0.73–0.95)	6.10E-03	0.76 (0.60–0.96)	0.019
			SBCGS	0.94 (0.86–1.03)	0.21	0.98 (0.83–1.16)	0.81
			meta	0.90 (0.85–0.95)	2.1E-04	0.84 (0.78–0.91)	7.3E-06
rs11546155	33,451,148	A/G	BPC3			1.19 (1.04–1.35)	9.1E-03
	159,054		CGEMS	1.12 (0.94–1.34)	0.21
			AABC	1.14 (1.02–1.27)	0.023	1.32 (1.14–1.54)	3.3E-04
			ROOT	1.11 (0.96–1.28)	0.15	1.18 (0.93–1.51)	0.18
			SBCGS	1.16 (0.98–1.38)	0.089	1.26 (0.94–1.70)	0.13
			meta	1.13 (1.05–1.21)	7.0E-04	1.23 (1.13–1.35)	2.0E-06

* NCBI 37 and from transcription starting site of TP53INP2

† rs8116198 is monomorphic in Asian population.

None of the tests for heterogeneity across studies was significant.

OR, odds ratio; CI, confidence intervals; ER, estrogen receptor

Table 4

Regulatory element annotation of variants that predicted expression of TP53INP2 using HaploReg [33].

Variant	Position*	Promoterhistone marks†	Enhancerhistone marks†	DNAse hypersensitivity	Proteinsbound	Motifschanged by the variant
rs1205339	-367,127		6 tissues including breast and blood (HMEC, MYO, HMEC.35)			ATF2, Mef2, Pax-4, Pou1f1, TATA
rs4911154	-295,993		Liver	Blood, liver	TCF4	RFX5
rs8116198	-177,893		24 tissues including breast and blood (HMEC, MYO, HMEC.35)		POL2, TBP, TR4	Rad21
rs6058107	-3,548	24 tissues including breast and blood (HMEC, MYO, HMEC.35)‡		28 tissues including breast and blood‡		AP-1,NF-E2
rs6060047	75,306		Multiple tissue types including blood and breast (HMEC, MYO, HMEC.35)‡	Multiple tissue types including blood		BATF, GCNF, HNF1, Irf, STAT
rs11546155	159,054	2 tissue types	4 tissue types			NRSF, Pou5f1, RXRA, Sin3Ak-20

* base pair from the transcription starting site of TP53INP2

† Normal mammary or breast cancer cell lines are indicated in parenthesis. HMEC.35: breast variant human mammary epithelial cells; MYO: breast myoepithelial primary cells; HMEC: mammary epithelial primary cells (vMHEC)

‡ Variants in strong linkage disequilibrium

Table 5

GAME-ON replication for SNPs related to the TP53INP2 gene.

	Test/ref allele	Study phase	Overall		ER-negative*
SNP			OR (95% CI)	P	OR (95% CI)	P
rs11546155	A/G	U4C	1.13 (1.05–1.21)	7.0E-04	1.28 (1.14–1.44)	5.0E-05
		GAME-ON	1.05 (0.99–1.10)	0.11	1.13 (1.03–1.25)	0.013
rs1205339	G/A	U4C	1.11 (1.05–1.18)	4.1E-04	1.22 (1.11–1.35)	8.0E-05
		GAME-ON	1.02 (0.98–1.07)	0.3	1.09 (1.01–1.17)	0.021
rs4911154	A/G	U4C	1.13 (1.06–1.20)	1.6E-04	1.25 (1.12–1.39)	5.5E-05
		GAME-ON	1.02 (0.98–1.07)	0.31	1.09 (1.01–1.17)	0.02
rs6058107	C/T	U4C	0.90 (0.85–0.95)	1.5E-04	0.90 (0.82–0.98)	0.021
		GAME-ON	0.96 (0.92–1.00)	0.045	0.91 (0.85–0.97)	5.0E-03
rs6060047	G/T	U4C	0.90 (0.85–0.95)	2.1E-04	0.82 (0.75–0.91)	1.2E-04
		GAME-ON	0.96 (0.91–1.00)	0.066	0.90 (0.84–0.97)	7.7E-03
rs8116198	G/A	U4C	0.93 (0.81–1.07)	0.33	1.04 (0.78–1.40)	0.79
		GAME-ON	0.97 (0.91–1.03)	0.28	0.94 (0.84–1.05)	0.28

*The overlapping study (BPC3) was removed from the meta-analysis in the discovery phase (U4C).

OR, odds ratio; CI, confidence interval; ER, estrogen receptor

Table 6

dbGaP datasets used in the our gene level expression-based GWAS analysis.

Accession Number	Study Name	Acronym	Breast Cancer Phenotype	Population
phs000851	African American Breast Cancer GWAS	AABC	3016 cases, 988 ER- cases, 2745 controls	African American
phs000812	Breast and Prostate Cancer Cohort Consortium GWAS	BPC3	1998 ER- cases, 3263 controls	European American
phs000147	Cancer Genetic Markers of Susceptibility Breast Cancer GWAS	CGEMS	1142 cases, 1145 controls	European American
phs000383	GWAS of Breast Cancer in the African Diaspora	ROOT	1657 cases, 403 ER- cases, 2029 controls	African American, African, African Barbadian
phs000799	Shanghai Breast Cancer Genetic Study	SBCGS	2790 cases, 490 ER- cases, and 2176 controls	Asian (Chinese)