Integrative Genomic Analysis Predicts Causative Cis-Regulatory Mechanisms of the Breast Cancer–Associated Genetic Variant rs4415084

Genome-wide association studies (GWAS) have identified thousands of common genetic variants associated with various traits and diseases, including cancer (1). However, most of these variants lie in noncoding regions of the genome where a direct link to gene function or regulation is difficult to assess (2). Furthermore, it is often unclear whether the true molecular perturbations associated with carcinogenesis, cancer progression, or therapeutic response indeed lie in the reported GWAS variants themselves or some other linked genetic variants. As a result, discovering the direct functional consequences of genetic variation at GWAS loci has been a critical missing step in utilizing the rich GWAS results to advance cancer research. We here address this important challenge by presenting an integrative computational framework that can facilitate the rapid identification of candidate functional regulatory variants in open chromatin regions and their target genes.

Some approaches are currently known for investigating candidate target genes and causative SNPs of GWAS variants, but they usually yield many false positives (3, 4). For identifying candidate target genes, two popular approaches correlate gene transcription level with the variants across a population of patients: one is expression quantitative trait loci (eQTL) analysis and the second is allele-specific expression (ASE) analysis. Both methods make use of genotype and mRNA transcription profiles in large patient cohorts available from various databases such as The Cancer Genome Atlas (TCGA). The first method has been successful for identifying genes that are correlated in overall mRNA level with the GWAS variant genotypes. For example, Li and colleagues (5) performed eQTL analysis on estrogen receptor–positive (ER⁺) breast cancer data and found SNPs correlating with the expression of essential genes, such as ESR1 and c-MYC. However, traditional eQTL analysis is only correlative and cannot distinguish between direct and secondary target genes. The ASE method uses patients who are heterozygous at a given GWAS SNP and tests for direct cis-regulated genes showing an imbalance in transcription activity between the two separate chromosome copies harboring different alleles of the variant. Such an imbalance is typically assessed from RNA sequencing (RNA-seq) data by counting the allele-specific coverage of exonic SNPs present in candidate genes within a certain distance from the GWAS SNP. Most studies to date, however, use unphased ASE analysis, that is, the exonic SNPs showing an allelic skew are not phased with the GWAS variant (5–7), thereby losing key information about whether the GWAS variant attenuates or promotes target gene transcription. In limited studies utilizing phased ASE (8), phasing accuracy has not been rigorously evaluated; as a result, testing for the consistency of phased ASE across patients may lose power when incorrectly phased haplotypes are used for some patients.

Our integrative approach improves upon these ideas in several ways. First, we restrict eQTL analysis only to those genes that lie within the same topologically associated domain (TAD) as a given GWAS SNP. Because TADs are thought to represent physical chromatin loops containing regulatory elements and their target genes (9), restricting the analysis to GWAS locus–containing TADs should help remove false positive genes and improve upon the current practice of choosing an arbitrary distance cutoff (5). We then apply ASE analysis on exonic SNPs properly phased with GWAS SNPs by removing patient samples whose reconstructed haplotypes show unstable phasing in our simulation study. This filtering method, which we termed the local-chromosome allele-specific expression (LCASE), effectively analyzes ASE on a local chromosome segment that is in LD with a GWAS SNP and increases statistical power by removing incorrect haplotypes from consideration. To find the causal regulatory SNP linked to a given GWAS SNP, we search for linked SNPs residing in functional open chromatin regions and computationally assess the SNPs' effect on candidate transcription factor (TF)–binding motifs, as we have previously done to show that the two well-known point mutations in the TERT promoter create new binding sites of GABP to reactivate TERT transcription (10).

We demonstrate the utility of our method by applying it to a breast cancer susceptibility region in 5p12, which is a GWAS hotspot harboring three noncoding GWAS SNPs replicated in previous studies (11–13). First, we show that all three GWAS SNPs may be targeting the same genes, the protein-coding gene MRPS30 and lncRNA RP11-53O19.1, both of which have been implicated in cancer (12, 14, 15). MRPS30 encodes a member of the mitochondrial ribosomal large subunits (14), suggesting the risk SNP's role in modulating mitochondrial activities. The lncRNA RP11-53O19.1, also known as breast cancer–associated transcript 54 (BRCAT54), is overexpressed in luminal A breast cancer subtype (ER⁺; ref. 15), suggesting its specific role in ER⁺ breast cancers. We then propose that an intergenic SNP, in LD with one of the 5p12 GWAS SNPs, is the predicted functional SNP. We provide multiple lines of evidence supporting that the risk allele of the predicted functional SNP increases the binding affinity of GATA3, an important TF known to cooperate with ESR1 and FOXA1 in ER⁺ breast cancers (16). Although this article focuses on the ER⁺ breast cancer risk hotspot in 5p12, the described method for dissecting the functional consequences of GWAS variants can be generalized to other noncoding loci implicated in a variety of cancers or traits.

Germline and tumor genotypes at tag SNPs from 979 patients were downloaded from the TCGA Data Portal and subsequently used for imputation. The following data were obtained from the NCI Genomic Data Commons (GDC) Legacy Archive (17): 693 patients' tumor copy number segmentation data, 788 patients' processed tumor gene expression data, and 795 ER⁺ breast cancer patients' tumor RNA-seq raw reads. Data from normal tissue and primary tumor were matched on the basis of patient barcodes. Patients having all four types of data available were selected to obtain a final set of 679 female ER⁺ breast cancer patients.

TCGA germline genotypes of tag SNPs measured by Affymetrix human SNP array 6.0 were mapped to 1000 Genomes Project phase 3 (18) variants. For each patient, genotypes with low quality (5th quantile in confidence score) were excluded. Because none of the three GWAS SNPs (rs4415084, rs10941679, rs7116600) were directly genotyped by the SNP array, imputation was performed in the 5p12 region using IMPUTE2 (19, 20) and the 1000 Genomes Project Phase 3 (October 2014 version; ref. 18) as a reference panel. Imputed genotypes were retained if the maximum genotype probability exceeded the threshold 0.9 and the minor allele frequency (MAF) exceeded 0.01 (for high imputation quality; ref. 20). To obtain a more accurate picture of germline genotypes, the genotypes from normal blood rather than primary tumor were used in all analyses to avoid miscalls from potential genotyping errors and somatic mutations.

Hi-C data were obtained from GSE66733 (21) and GSE51687 (22) for MCF-7, and from GSE53463 (23) and the Encyclopedia of DNA Elements (ENCODE; ref. 24) ENCSR549MGQ for T-47D. MCF-7 TAD structures were determined on the basis of insulation score from Hi-C data (GSE66733). The TAD information in T-47D was obtained from the ENCODE dataset ENCFF075QYD. We used the MCF-7 TAD information in the analysis, as the two datasets were similar. MCF-7 ChIA-PET interactions mediated by CTCF, ESR1, and POLR2A were obtained from ENCODE (24).

We performed TCGA-based eQTL analysis by constructing a multivariate linear model that regressed the expression level of each gene against the copy number of the gene and the genotype status at a given GWAS SNP locus. For the ith GWAS SNP in 5p12, we encoded its germline genotype G_i as the number of risk alleles (G_i ∈{0, 1, 2}). To determine the copy number of the |j$|th gene residing within the 5p12 TAD region, we took length-weighted average (⁠|\overline {{s_j}} $|⁠) of tumor copy number segment data covering the gene and transformed it into copy number |C{N_j} = 2*{2^{\overline {{s_j}} }}$|⁠. The RNA-seq RPKM value of gene j was log-transformed as |{E_j} = lo{g_2}( {RPK{M_j} + 1} )$|⁠. Multivariate linear regression (25), |{E_j}_\ = {\alpha _{ij}} + {\beta _{ij}}\ {G_i}{\bm{\ }} + {\gamma _{ij\ }}C{N_j} + {\epsilon _{ij}}$|⁠, was then performed for each pair of |i \in \{ {GWAS\ SNPs\ in\ 5p12} \}$|and |j \in \{ {genes\ in\ 5p12\ TAD} \}$|⁠, where |{\alpha _{ij}}$| is the intercept and |{\epsilon _{ij}}$| the error term. Genes with |\overline {RPKM} \ge 1$| (⁠|\overline {RPKM} $|⁠: mean expression among tumor samples) and |p_{ij}^G \le \alpha $| (⁠|p_{ij}^G:\ p$|-value of |{\beta _{ij}}\ $|⁠; α = 0.05) were called as significant eQTL target genes. Bonferroni correction for multiple hypothesis testing was further applied using |p_{ij}^G \le {\frac{\alpha }{n}$|⁠, where n is the total number of genes tested in the TAD (n = 22, thus |{\frac{\alpha }{n}} = 0.0023$|⁠).

We tested for within-patient differential transcription between the two chromosome copies harboring different alleles of a given GWAS SNP using a modified version of ASE analysis, which we termed LCASE. To distinguish between the transcripts from different chromosome copies, we utilized exonic SNPs (GENCODE v19; ref. 26) that were heterozygous based on imputation results. For each exonic SNP, we selected a subset of patients with heterozygous genotypes at both the GWAS SNP and the exonic SNP. For each patient in this subset, we used SHAPEIT2 (27) to determine the phase between the GWAS SNP and exonic SNP. We denoted the risk and protective alleles at the GWAS SNP as risk and pro, respectively, and denoted the reference and alternative alleles of the exonic SNP under consideration as 0 and 1, respectively. To remove uncertain phasing, we further sampled haplotypes n times (n = 100) from the haplotype model generated by SHAPEIT2 and computed the occurrence probability |( {{p_{risk|0}},\ \ {p_{risk|1}}} )$| of the two phases, thereby removing patient samples with no dominant phase as indicated by the condition |\max ( {{p_{risk|0}},\ \ {p_{risk|1}}} ) \ \lt\ 0.9$| (Supplementary Fig. S1). We then counted the allelic-specific coverage of each exonic SNP from high-quality RNA-seq reads (MAPQ ≥ 20) and tested for transcription imbalance using two statistical tests depending on the sample size K: binomial test (when K < 5) and Wilcoxon signed-rank test (when K ≥ 5).

A list of GWAS variants in the 5p12 region was obtained from the NHGRI GWAS Catalog (28), and the risk alleles were extracted from Stacey and colleagues (11) and Thomas and colleagues (13). All common (MAF ≥ 0.01) SNPs from 1000 Genomes Project Phase 3 in high LD (⁠|{r^2} \ge 0.8$|⁠) with any of the 5p12 GWAS SNPs were selected for prioritization. All r² values were calculated on the basis of the 1000 Genomes Project Phase 3 EUR population. For SNP prioritization, DNase I hypersensitivity sites (DHS) in MCF-7 and T-47D were collected from the ENCODE (24) database, including DHS under estradiol treatment (Supplementary Table S1). Only those LD SNPs that overlapped with a narrow peak from at least one of DHS replicate experiments were considered for further investigation. We also used histone modification and TF-binding profiles to facilitate functional SNP prioritization. Histone modification data were obtained from ENCODE and the Gene Expression Omnibus (GEO) database (GSE26831, GSE63109, and GSE69112). Chromatin immunoprecipitation sequencing (ChIP-seq) data in MCF-7 and T-47D for TFs, including ESR1, PGR, GATA3, and FOXA1 known to be important in breast cancer, were obtained from ENCODE and GEO (Supplementary Table S2). Raw ChIP-seq data from GEO were mapped using BWA (default parameters; ref. 29), and peaks were called using MACS2 (TF ChIP-seq: default parameters; histone modification ChIP-seq: broadpeak mode, 0.1 FDR; ref. 30).

We collected TF position weight matrices (PWM) from HOCOMOCO Human v10 (31), FACTORBOOK (32, 33), TRANSFAC (34), JASPAR vertebrates (35), and Jolma2013 (Supplementary Table S3; ref. 36). To identify motifs potentially affected by a SNP, we used the program FIMO (version 4.12.0; ref. 37) to scan the two 50-bp sequences centered at the different alleles of each candidate SNP (FIMO threshold |{10^{ - 3}}$|⁠). We then sought for candidate SNP-TF pairs for which the sequence harboring one allele of the SNP matched the TF motif, whereas the sequence harboring the other allele did not. For each SNP–TF pair, the statistical significance of creating or disrupting the motif via the SNP was evaluated with two approaches: mutation simulations and, when available, ChIP-seq k-mer enrichment analysis. The first approach compared the change in motif score caused by the SNP against the changes caused by random single-nucleotide mutations simulated 5,000 times, where the motif score was defined as |\mathop \sum_{l = 1}^L \log {p_{l,{n_l}}}\ $|(⁠|{p_{l,{n_l}}}:$| the probability of nucleotide n_l at position l in the PWM of length L; Supplementary Methods). The second approach tested for the differential enrichment of k-mers (e.g., k = L) harboring different alleles by comparing their occurrence frequency in TF ChIP-seq peaks against that in control regions (random open chromatin regions in MCF-7 or T-47D; χ² test; Supplementary Methods; ref. 38).

The candidate SNP–TF pairs from motif analysis were further filtered on the basis of the correlation structure between each candidate TF and its eQTL target gene. A list of human TFs was obtained from AnimalTFDB (39), and only the expressed TFs (⁠|\overline {RPKM} \ge 1$| among ER⁺ breast cancer patients) overlapping with the candidates from motif analysis were further considered and divided into three groups: candidate activators (Pearson correlation coefficient |r \ge {r_{ + 2\sigma }}$|⁠), candidate repressors (⁠|r \le {r_{ - 2\sigma }}$|⁠), and uncorrelated TFs (⁠|{r_{ - 2\sigma }}\ \lt\ r\ \lt\ {r_{ + 2\sigma }}$|⁠), where the thresholds |{r_{ - 2\sigma }}$| and |{r_{ + 2\sigma }}$| were chosen based on the distribution of correlation coefficients between random pairs of TF and non-TF genes (⁠|{r_{ - 2\sigma }} = - 0.337$|⁠; |{r_{ + 2\sigma }} = 0.409$|⁠; Supplementary Methods). We then narrowed the list of candidate TFs based on four models shown in Supplementary Fig. S2. In this article, we focused on an example of enhancer-activating risk allele. In this case, where the target gene expression was positively correlated with the number of GWAS risk allele, we selected TFs passing the following filters as candidate activators: (i) the risk allele of the candidate enhancer SNP significantly increased the motif score of the TF; (ii) the TF expression level was positively correlated (⁠|r \ge {r_{ + 2\sigma }}$|⁠) with the target gene expression level among patients having two copies of the risk allele; (iii) the TF-target gene expression correlation in patients carrying two copies of the risk allele (two copies of the motif) was stronger than that in patients carrying no risk allele (no motif).

For TFs with ChIP-seq data available in MCF-7 and T-47D, we tested for allele-specific binding (ASB) by searching for a skew in the read coverage of candidate heterozygous SNPs located within peaks. The significance of a skew was measured by using the binomial test assuming no bias. To determine the genotype status, tag SNP data for MCF-7 and T-47D were obtained from CCLE project (40), and imputation was performed using the same pipeline described above. To check that the candidate loci have no deletion or amplification events, which might complicate the ASB analysis, copy number data for the cell lines were obtained from CCLE (40) and ENCODE (GSE40698; ref. 24). ChIP-seq reads from different replicates were combined and deduplicated.

We integrated multiple genomic analyses to identify systematically the triplets of functional SNP, corresponding TF regulator, and target gene (Fig. 1A and B). Our main hypothesis was that noncoding SNPs modulating cancer risk would likely lie inside putative regulatory elements to promote or attenuate the binding affinity of relevant TFs, which would be reflected at the mRNA level of the target gene. A schematic representation of our framework is shown in Fig. 1A. Target genes were identified by using within-TAD eQTL and LCASE analyses of GWAS SNPs, and functional SNP-TF candidates were determined by integrating motif analysis, expression correlation analysis, and epigenetic annotation of regulatory regions (Materials and Methods) for GWAS SNPs and their high LD SNPs.

We applied our computational framework to a breast cancer susceptibility region in 5p12. As there are three ER⁺ breast cancer GWAS SNPs in the 5p12 risk hotspot (Fig. 2; refs. 11, 13), we first performed eQTL analysis for individual SNPs to investigate whether they all targeted the same genes. We restricted the analysis to the TAD containing all three GWAS SNPs based on the MCF-7 TAD information (21). The 5p12 TAD region (chr5:43,480,001-45,000,000, hg19) containing all three GWAS SNPs encompassed 7 protein-coding genes and 15 noncoding RNAs (GENCODE v19). Among them, the protein-coding gene MRPS30 and lncRNAs RP11-53O19.1 and RP11-53O19.3 showed significant eQTL correlation (genotype effect P ≤ 0.05) with all three GWAS SNPs, but RP11-53O19.3 had only borderline significance with rs7716600 (Supplementary Tables S4 and S5). Notably, the GWAS SNP rs4415084 correlated strongest with these genes compared with other GWAS SNPs (MRPS30, |P = 1.39 \times {10^{ - 5}}$|⁠; RP11-53O19.1, |P = 6.19 \times {10^{ - 6}}$|⁠; RP11-53O10.3, |P = 1.23 \times {10^{ - 4}};\ $|eQTL analysis in Materials and Methods; Fig. 3A). Furthermore, rs4415084 was the only GWAS SNP that remained significant after Bonferroni correction (Supplementary Table S4).

It should be noted that the three genes are likely coexpressed, as RP11-53O19.1 and MRPS30 share a divergent promoter and RP11-53O19.3 is contained in several elongated transcripts of MRPS30 (Supplementary Fig. S3). As RP11-53O19.3 might be a processed transcript of MRPS30, our subsequent analyses focused only on MRPS30 and RP11-53O19.1. The positive correlation between the number of risk alleles and the expression level of these target genes was also confirmed in both isoform-level eQTL analysis in TCGA data (Supplementary Table S6) and gene-level eQTL analysis in GTEx normal breast tissue (MRPS30, |P = 5.46 \times {10^{ - 3}}$|⁠; RP11-53O19.1, |P = 4.14 \times {10^{ - 4}}$|⁠; Supplementary Fig. S4). Together, these results showed that all three GWAS SNPs in the 5p12 region may target the same set of genes within the TAD, with rs4415084 displaying the strongest effect.

Because eQTL analysis is only correlative and cannot distinguish between trans-regulation and cis-regulation, we sought evidence for cis-regulation by utilizing heterozygous exonic/UTR SNPs properly phased with the GWAS SNPs (Fig. 3B; Materials and Methods). We first examined all common exonic variants (MAF ≥ 0.01, 1000 Genomes Phase 3, EUR) in the protein-coding gene MRPS30 and found three SNPs suitable for LCASE analysis (rs61754779 and rs34522103 in exons; rs79210252 in 3′UTR; Fig. 3C). There were 5 patients with heterozygous genotypes at both the GWAS SNP rs4415084 (C/T) and the SNP rs61754779 (G/C) located in the first exon. Among these five, two were removed because of uncertain phasing. The risk haplotypes in these 3 patients were all rs61754779-G|rs4415084-T; equivalently, all of the protective haplotypes were rs61754779-C|rs4415084-C. RNA-seq read coverage at this exonic SNP displayed a significantly higher transcription activity from the risk haplotype in 2 patients (⁠|P = 1.33 \times {10^{ - 4}}$|⁠, |P = 9.72 \times {10^{ - 17}}$|⁠; one-sided binomial test; Fig. 3C), but the skew in the third patient was not significant (P = 0.64). Similar analyses at the exonic SNPs rs34522103 and rs79210252 identified a different set of three patient samples that showed the same direction of transcription imbalance toward the risk haplotype, and the skew was statistically significant in two of these patients (rs34522103: |\ P = 6.71 \times {10^{ - 47}}$|⁠; rs79210252: |P = 1.22 \times {10^{ - 3}}$|⁠, P = 0.28; one-sided binomial test; Fig. 3C).

Although the LCASE in MRPS30 was significant in several patients, the sample size was small, because exonic SNPs were generally rare in the genome. We thus sought for more evidence in MRPS30 ASE by analyzing additional SNPs in the noncoding RNA (ncRNA) located at the 3′ end of MRPS30 (GENCODE track, Supplementary Fig. S3). De novo assembly of transcripts using StringTie (41) on tumor RNA-seq reads revealed multiple MRPS30 elongated transcripts covering both MRPS30 3′ ncRNA and RP11-53O19.3 (Supplementary Fig. S3), suggesting that these two ncRNAs could arise from abnormal elongation of the protein-coding MRPS30 gene. This elongation hypothesis was also supported by high correlation in expression level between MRPS30 and RP11-53O19.3 (r = 0.95).

We therefore tested whether there were additional LCASE SNPs falling in this elongated region and found four transcribed SNPs that could be phased accurately in a large number of patient samples (n = 23, n = 248, n = 238, n = 240, respectively). Consistent with the analysis of exonic SNPs, all of these SNPs in the elongated region displayed significantly higher transcription activity on the risk chromosome (⁠|{P = 9.5 \times {10^{ - 7}},\ P = 5.5 \times {10^{ - 13}},\ P = 7.7 \times {10^{ - 17}},\ P = 3.2 \times {10^{ - 14}},$| respectively; Wilcoxon signed-rank test; Fig. 3D; Supplementary Fig. S5). Similar analyses in RP11-53O19.3 found that five out of eight transcribed SNPs showed a skew toward the risk chromosome, while two other SNPs close to the 3′ end displayed an opposite trend (Supplementary Fig. S6). Together, the eQTL and LCASE analyses thus demonstrated that the breast cancer risk allele rs4415084-T was associated with increased MRPS30 and RP11-43O19.1/3 expression among ER⁺ breast cancer patients in a cis-regulating manner (Fig. 3C and D; Supplementary Figs. S5–S7).

We next prioritized candidate LD SNPs using epigenetic information. We gathered 123 variants in high LD (⁠|{r^2} \ge 0.8$|⁠, 1000 Genomes Phase 3, EUR) with any of the three GWAS SNPs rs4415084, rs10941679 or rs7716600, including the GWAS SNPs themselves (Fig. 2). Using the DNase I hypersensitivity (DHS) data from ENCODE in two breast cancer cell lines (MCF-7 and T-47D; Fig. 2; refs. 24, 42), we found that 16 of these SNPs resided in DHS peaks, all of them being in high LD with rs4415084. Among the 16 SNPs, three SNPs were located near the center of DHS peaks (Supplementary Fig. S8) and were thus prioritized for further investigation in motif analysis and TF correlation analysis.

Analyzing the prioritized candidate SNPs provided multiple lines of evidence for rs4321755 [chr5:44646195; about 163 kb upstream of MRPS30 transcription start site (TSS)] being the functional SNP. The risk allele of SNP rs4321755 could be determined to be T, because this SNP was in high LD with GWAS rs4415084 (⁠|{r^2} = 0.988,\ $|1000 Genomes Phase 3, EUR), and rs4321755-T was the allele in phase with the GWAS risk allele rs4415084-T. First, our motif analysis showed that the risk allele rs4321755-T created a GATA3-binding motif; the 9 bp sequence covering rs4321755 matched a GATA3 motif only when it harbored the allele rs4321755-T, and not when it contained the alternative allele rs4321755-C (Fig. 4A, FIMO |P = 1.59\ \times {10^{ - 4}}$|⁠). We further confirmed that the T-to-C conversion at the SNP resulted in a dramatic frequency drop from 93.7% to 0.0% in the GATA3 PWM, which was a significant change based on both mutation simulations (P = 0.003; Supplementary Methods), and ChIP-seq k-mer enrichment analysis (P = 0.025, χ² test; Supplementary Methods).

Second, GATA3 showed a consistent pattern in the TF-target gene correlation analysis: among the TCGA patients, GATA3 expression positively correlated with the expression levels of predicted target genes MRPS30, RP11-53O19.1, and RP11-53O19.3, supporting GATA3′s role as a transcription activator (Fig. 4B). Furthermore, the positive correlation was strongest (and above |{r_{ + 2\sigma }}$| for MRPS30 and RP11-53O19.3) in patients carrying the rs4321755-T/T genotype, moderate in patients with the rs4321755-C/T genotype, and weakest in patients with the rs4321755-C/C genotype (Fig. 4B), consistent with our model of enhancer-activating risk allele (Supplementary Fig. S2). This correlation trend implied that the regulation mediated by GATA3 was strongest when two copies of GATA3 motif existed at the SNP on both paternal and maternal chromosomes, while the regulation was weakest when both GATA3-binding sites at the SNP were disrupted.

We next sought direct experimental evidence of GATA3 binding at the predicted causal SNP. Consistent with our prediction, we found that rs4321755 was located at the center of a peak in the T-47D GATA3 ChIP-seq data from ENCODE (q-value |= {10^{ - 4}},$| 2bp from the summit; Fig. 4C). Moreover, FOXA1 and PGR ChIP-seq in T-47D also showed cobinding with GATA3 near rs4321755, supporting that this SNP indeed lies within an active regulatory region in breast cancer, given that FOXA1 is a pioneer factor in ER⁺ breast cancers (43) and PGR can interact with ER (Fig. 4C; ref. 44); ESR1 ChIP-seq showed no direct binding of ER itself near rs4321755 (GSE32465; ref. 24). Furthermore, this regulatory region was evolutionarily conserved (100 vertebrates base-wise conservation by PhyloP; Fig. 4D; ref. 45), indicating its important role retained through evolution. To measure GATA3 ASB, we imputed the genotypes and found T-47D to be heterozygous at rs4321755. The copy number of the segment containing the SNP suggested that T-47D had no deletion or amplification at this locus, informing that the null hypothesis should be unbiased binding between the two alleles (Materials and Methods). We obtained and deduplicated DNase-seq and multiple TF ChIP-seq reads for ASB analysis in T-47D. As shown in Fig. 4E, GATA3 ChIP-seq reads contained more of the risk allele rs4321755-T than the protective allele rs4321755-C (P = 0.019, one-sided binomial test). DNase-seq and PGR ChIP-seq also showed significantly more reads with T than C (DNase-seq, P = 0.032; PGR, P = 0.002; one-sided binomial test); FOXA1 had a similar imbalance, although the P value was not significant due to low coverage at the SNP (Fig. 4C and E). Quite interestingly, MCF-7, which was homozygous for the protective allele C at rs4321755, had a closed chromatin configuration at the SNP (Supplementary Fig. S8). Thus, the SNP likely contributed to the lack of GATA3 binding in MCF-7, potentially jointly with GATA3′s impaired DNA-binding ability caused by a heterozygous frameshift mutation in MCF-7 (46). Together, these findings provided experimental evidence that the putative causal SNP rs4321755 might not only influence GATA3 ASB, but also modulate differential chromosomal accessibility and, thus, enhancer activity.

We next explored available chromatin interaction data to assess potential long-distance enhancer–promoter interactions, as the candidate SNP rs4321755 was about 163 kb away from the MRPS30 TSS. However, rs4321755 was in DHS only in the T-47D cell line, for which no GATA3 or other ChIA-PET data were available. In MCF-7, data were available for Hi-C (21) and ChIA-PET of breast cancer–related factors, such as ESR1 and CTCF (47). However, no significant chromatin interactions originating from the regulatory region were observed in MCF-7, probably because this enhancer was specific to T-47D and was not accessible in MCF-7, which carried the C/C genotype at rs4321755 (Supplementary Fig. S8). Comparing the Hi-C data in MCF-7 versus T-47D at 40 kb resolution using edgeR (48), we found a significantly higher contact count in the heterozygous T-47D than the homozygous MCF-7 between the bin containing rs4321755 and the bins spanning MRPS30 (edgeR |FDR\ \lt\ 1 \times {10^{ - 2}}$|⁠) and RP11-53O19.1 (edgeR |FDR\ \lt\ 1 \times {10^{ - 4}}$|⁠; Supplementary Methods; Supplementary Table S7).

This study presented a computational framework for systematically investigating the functional consequences of GWAS SNPs. We applied our method to a breast cancer susceptibility region in 5p12 and discovered a causative SNP with multiple lines of evidence supporting its function in modulating GATA3-binding affinity. This causative SNP potentially explained the molecular mechanism of one 5p12 GWAS SNP, rs4415084, while the other two GWAS SNPs in the 5p12 region had no candidate SNPs found in open chromatin. It is currently unclear whether these three GWAS SNPs function through the same or independent regulatory elements. For instance, the risk allele G of rs10941679 was completely contained in the background of the risk allele rs4415084-T that had higher allele frequency and thus broader impact in population (233/234 patients with rs10941679-G carry rs4415084-T; rs10941679 MAF: 0.23; rs4415084 MAF: 0.41; 1000 Genomes Project Phase 3, EUR; ref. 18). Moreover, Stacey and colleagues (11) showed that the risk in 5p12 could be explained by either rs4415084 or rs10941679, with the significance of rs4415084-T remaining after correcting for rs10941679-G and vice versa. Although a previous study found that the GWAS SNP rs10941679 itself could be causal (12), this particular SNP was not in open chromatin regions of the breast cancer cell lines that we examined (Supplementary Fig. S8). Another study reported correlation between our target gene MRPS30 and rs7716600 genotype, but no candidate causative SNPs in TF-binding sites were discussed (49). Although the breast cancer susceptibility harbored in 5p12 is not totally understood, we here propose one functional SNP that may directly link rs4415084 to the regulation of predicted target genes.

Our integrative approach facilitated the identification of putative target genes by combining eQTL and LCASE analysis. We found breast cancer–specific TADs containing GWAS SNPs and performed intra-TAD eQTL analysis. As TADs are thought to provide physical subdivisions of cis-regulation, this approach should help reduce the number of false positives that currently challenge the traditional eQTL analysis. In our LCASE analysis, we controlled for phasing quality by removing samples that could not be phased confidently upon simulation, enabling a direct measurement of differential transcription activity from two local chromosome copies harboring different GWAS alleles. The target genes identified by our method were already implicated in cancer: the protein-coding gene MRPS30, closely related to mitochondrial activity, was also linked to the 5p12 risk variants (rs10941679 and rs7716600) in earlier studies (12, 49). Similarly, the candidate target lncRNA RP11-53O19.1 was previously found to have significantly higher expression in luminal A (ER⁺) breast cancers compared with other subtypes (12). Even though the expression levels of these target genes may be directly modulated by the identified SNP, cancer is a complex multigenic disease, and understanding how the altered expression levels contribute to cancer predisposition requires further investigation, likely involving a systems-level approach.

Our framework predicted candidate (causative SNP, TF, target gene) triplets to address the challenging problem of discovering functional SNPs that, we hypothesized, might disrupt TF-binding activities (Fig. 5). We integrated cell type–specific epigenetics profiles, motif analysis, and expression correlation signatures to reduce false positives. Our method actually yielded a ranked list of candidate TFs. Our top candidate GATA3 was validated through GATA3 ChIP-seq data and ASB analysis. Besides GATA3, TCF7L1, NR3C1, and ETS1 also qualified as candidate TFs when the TF-gene correlation thresholds were chosen to be less stringent (Supplementary Fig. S8). These predictions would benefit from future ChIP-seq data as they become available and also from imputing binding profiles (50). Currently, a key challenge in annotating functional regulatory elements is the lack of 3D chromatin interaction data that can help validate the physical interaction of predicted enhancer–promoter pair. For this study, GATA3 ChIA-PET data in T-47D, currently missing, would be most suitable for validating the predicted interaction between the identified enhancer and MRPS30 promoter.

We present an integrative computational method using genomic and epigenomic data to identify causative regulatory variants that may directly modulate cancer predisposition. Application of our method to a breast cancer susceptibility region in 5p12 reveals the intergenic SNP rs4321755, in LD with the GWAS SNP rs4415084, as the candidate causative variant. We propose that the risk allele rs4321755-T significantly increases GATA3-binding affinity and therefore results in upregulation of the predicted target genes MRPS30 and RP11-53O19.1. Our computational framework can be extended to investigate other genetic variants modulating cancer susceptibility, contributing to understanding new pathways in tumorigenesis and developing personalized prevention of cancer.

No potential conflicts of interest were disclosed.

The content is solely the responsibility of the authors and does not necessarily represent the official views of the NIH.

Conception and design: J.S. Song

Development of methodology: Y. Zhang, M. Manjunath, J.S. Song

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): Y. Zhang, M. Manjunath, J.S. Song

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): Y. Zhang, M. Manjunath, S. Zhang, D. Chasman, S. Roy, J.S. Song

Writing, review, and/or revision of the manuscript: Y. Zhang, M. Manjunath, S. Zhang, D. Chasman, S. Roy, J.S. Song

Study supervision: J.S. Song

The results appearing here are in part based upon the data generated by the TCGA Research Network (http://cancergenome.nih.gov/, dbGaP accession number phs000424.v6.p1 on 05/06/2016) and the GTEx Project, supported by the Common Fund of the Office of the Director of the NIH and by NCI, NHGRI, NHLBI, NIDA, NIMH, and NINDS. We acknowledge the ENCODE consortium that generated the datasets used in the article. Y. Zhang, M. Manjunath, and J.S. Song were supported by the grants 1U54GM114838, awarded by National Institute of General Medical Sciences (NIGMS) through funds provided by the trans-NIH (NIH) Big Data to Knowledge (BD2K) initiative, and R01CA163336. S. Zhang, D. Chasman, and S. Roy were supported by the NIH BD2K grant U54 AI117924.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.