Genetic variation in putative regulatory loci controlling gene expression in breast cancer

SNP Mining and Genotype Results.

The candidate genes for SNP analysis were selected from ≈4,000 MEDLINE entries and different databases (Online Mendelian Inheritance in Man and Human Genome Organization) to create the genotype profile of all genes involved in the regulation of the redox level in the cells, ROS-mediated signaling, and repair of DNA damage caused by ROS (Fig. 1). A total of 1,030 SNPs was selected by the SNP data mining approach described in Edvardsen et al. (9). Briefly, each of the candidate genes extracted from PubMed recourses was matched with the official Human Genome Organization gene name, and information about sequence variation was obtained. The annotation of each SNP was followed in several databases, and data on its frequency in the Caucasian population were obtained. Approximately 4,000 SNPs could be identified in the initial 233 candidate genes by computer annotation (Fig. 1A). SNP selection was performed based on putative gene function and SNP frequencies. We performed actual genotyping of 1,030 of these SNPs in 213 genes on chromosomes 1–X. Of these, 725 (≈69%) were successfully genotyped in the first round. One hundred four (≈14.3%) of the 725 genotyped SNPs were with frequency of <5%. An additional 38 SNPs were not in Hardy–Weinberg equilibrium (≈5%), leaving 583 SNPs for association studies. We studied the tumor genome-wide expression effect of these 583 SNPs that reside in 203 of the selected genes (1–19 SNPs per gene).

Statistical Analysis of SNP–mRNA Expression Associations.

Three complementary analyses [ANOVA, quantitative mutual information score (QMIS), and leave-one-out cross-validation (LOOCV)] (10–12) were used to analyze the data and to demonstrate possible parametric and nonparametric approaches to assess SNP–transcript association. We computed an association matrix whose entries represent the P values of association for every SNP/transcript pair (11). Let N and M denote the number of SNPs and transcripts, respectively. For each pair (s, t) of SNP and transcript, we computed an association score and a corresponding P value, P_st, using one of the three methods described (see Materials and Methods). The resulting N-by-M matrix P is called the association matrix. Given a significance threshold 0 < P < 1, we state that s is associated with t if P_st < P. SNPs whose corresponding rows have significantly many entries with low P values are potential regulators of expression levels. Subsets of SNPs and subsets of transcripts that belong to biclusters were identified by a method adapted from Tanay et al. (13) (Fig. 1B).

Significant Associations.

SNP expression associations with the best P values (P < 0.001) detected by both ANOVA and QMIS revealed regulatory SNPs in trans in 115 genes, such as AKT1, AKT2, CALM3, CDC25B, DPYD, FOS, IER3, IGF1R, IGF2, IGF2R, IL8, IL10, IL10RA, NFKB1, PRKCA, PPP3CA, PPP2R4, GCLM, TGFBR3, PPP1R2, and others (see the full list in Table 5, which is published as supporting information on the PNAS web site, including P values for all SNP–transcript pairs, where 10⁻⁶ was the lowest). For example, the presence of the variant G or C allele for SNPs in PPP1R2 and TGFBR3, respectively, led to an increased expression of QDPR and FANCA, even more pronounced in the homozygous variant genotype (Fig. 4, which is published as supporting information on the PNAS web site). Of the 115 genes, 43 genes harbored SNPs associated with mRNA expression both in cis, i.e., to the expression of their own gene, and in trans by both QMIS and ANOVA when requiring P < 0.001 in both methods (highlighted in Table 1; P values for each association are in Tables 5 and 6, which are published as supporting information on the PNAS web site). For example, SNPs that significantly associated with the expression of their own genes were found in GSTM3, XPC (two SNPs), GSTA4, GSTP1, ABCC1 (two SNPs), TYMS (four SNPs), CALM3, and MAPK1 (see the full list in Table 6). When we broadened the search for in cis interactions with mRNA coded 4 Mbp upstream and downstream of an SNP, 107 significant associations were found at P < 0.025 (Table 6; some exemplified in Fig. 5, which is published as supporting information on the PNAS web site). Some of the SNPs regulating in cis form are in strong linkage disequilibrium (LD).

For SNP–transcript pairs with strong associations supported on many samples we expected both QMIS and ANOVA P values to be small but not necessarily in the same range. For SNPs with small number of samples with a given genotype, we can see disagreement between the scores. ANOVA and QMIS P values were calculated based on different assumptions about the distributions of genotypes and expression values for each SNP association pair; these assumptions may not necessarily hold for a small number of samples. Therefore, restricting our analysis to SNP–transcript pairs significant with respect to both scores is an additional way to eliminate false positives. Table 2 shows the number of SNP–transcript pairs for various P value cutoffs and sizes of corresponding overlaps. The last column shows the expected overlap between randomly selected subsets of corresponding sizes. Note that the actual overlap in SNP–transcript pairs is much larger than was expected in a random case.

LOOCV analysis revealed cases in which the expression of only one transcript was sufficient to correctly classify the cases according to their genotypes. For example, the expression of TFF1 [trefoil factor 1 (breast cancer, estrogen-inducible sequence)] on 21q22.3 was sufficient to correctly predict all of the genotypes for an SNP (rs 2228001) in the XPC gene on 3p25, with a success rate of 0.9, and addition of other genes did not improve the classification. Similarly, the expression of YARS2 on 12p11.21 was sufficient to predict alone the grouping according to the genotype in rs881878 of EGF on 4q25 (Table 7, which is published as supporting information on the PNAS web site).

We identified SNPs in EGF, IL1A, MAPK8, XPC, SOD2, and ALOX12 genes whose associations to gene expression are significant under all three methods, namely ANOVA and QMIS with P <0.001 and classification success of LOOCV analysis >90%, indicating the presence of regulatory SNPs in these genes.

Statistical Overabundance of Significant Associations.

We applied FDR analysis (14) to compare the actual association matrix to random data. FDR measures the ratio of expected and observed numbers of SNP–transcript association pairs with a given score or better. The number of associations in the data were determined by counting the number of entries in the association matrix less than or equal to a given threshold. For QMIS, 769 SNP–transcript association pairs with P ≤ 1.0E-04 were observed. In random data one may expect to find only 150 such pairs, as inferred from QMIS P value and the size of the association matrix, which represents an FDR of 0.2 (Fig. 2). Note that one SNP can show significant association to more than one transcript. The graphical presentation of the distributions of observed and expected numbers of SNP–transcript pairs with a certain P value or lower for QMIS is shown in Fig. 2. Right shows the corresponding FDR. For the ANOVA scores (Fig. 2 Left), 571 SNP–transcript association pairs with P ≤ 1.0E-04 were observed, yielding a FDR of 0.6 as estimated from permuted data (Fig. 6A, which is published as supporting information on the PNAS web site) (15). The level of significance (P < 1.0E-04) may not seem to be stringent enough for this type of study with a huge burden of multiple tests. To address this issue and to assess the true significance of the results, we performed an overabundance analysis and estimated the distribution of ANOVA P values using simulations. To assess the baseline false-positive rate, 50 random SNPs in different genes were selected, and genotype data were permuted 100 times. For P ≤ 1.0E-05, FDR = 0.55, and for P ≤ 1.0E-04, FDR = 0.7. The observed and expected ANOVA P values are shown in Fig. 6B, and it can be seen that they are very similar to the results shown in Fig. 6A.

Master Regulators.

SNP loci with an exceptionally large set of significantly associated transcripts are putative regulators of many transcripts. They potentially affect the expression levels or the mode of operation of transcription factors or of noncoding RNA-mediated regulation, directly or indirectly. Thus, these SNPs affect the transcription or degradation rates and hence the expression levels of many transcripts. SNPs that had exceptionally dense rows in the association matrix, P, belong to the genes PRKCA, CALM3, CYP2C19, IGF1R, IGF2R, and XDH. An exact definition of row density is given in ref. 15. For some of the putative master-regulating genes (IGF2R and CALM3) associations both in cis and in trans were detected (Table 1).

Functional Grouping of Regulated Transcript, Pathway, and Gene Ontology (GO) Analysis.

The trans interactions observed here suggest both “positional” and “functional” explanations. The positional scenario may involve either LD or epigenetic events such as common domains of relaxed chromatin structure along the chromosomes (16). The functional scenario may involve interactions of the kind between receptors and their ligands, transcription factors, and genes under their control, the genes for which do not have to reside in vicinity to each other. Of that kind, we observed a strong association between SNPs in several growth factors and the expression of their receptors like EGF/EGFR, IL1/IL1R, and TGFB2/TGFBR1–3. In fact, in this study we used the genotypes from a given predefined pathway (the ROS metabolizing and signaling pathway) enriched for GO terms like ATP binding, phosphate metabolism, phosphorylation, and tyrosine kinase activity to search for associations from an unselected pool of whole-genome mRNA transcripts. If we inversely examine the GO terms of the mRNA that we find associated with these SNPs (Table 8, which is published as supporting information on the PNAS web site), we often find the same GO terms as for the candidate genes (Fig. 7, which is published as supporting information on the PNAS web site), suggesting that the observed associations are within the same functional pathway.

A Common Set of SNPs Associated with a Common Set of Transcripts (Biclusters).

A bicluster is a submatrix of the association matrix, that is, a subset of SNPs and transcripts. The goal of the bicluster analysis is to identify significantly dense biclusters, in which the participating SNPs share significantly many common associations with the participating transcripts (Fig. 1B, yellow rectangles). When several SNP loci, especially when not in LD, associate with the same transcript or group of transcripts, this may be evidence of this expression phenotype being a complex trait, affected by several genetic events. Several significant biclusters were found in the data, one of them consisting of four SNPs (in IL1B, IER3, and NOX3, all related to stress response) and 82 transcripts (see the full list in Table 8). The GO term “cytosolic small ribosomal subunit” was significantly overrepresented in this set of transcripts (Bonferroni corrected P = 0.002) together with other related GO terms (Fig. 1C). When the same set of loci commonly associate with a large number of transcripts, we have cross-confirmation of the individual associations as well as a possible multilocus effect on a pathway or a biological process.

Common Blocks of Regulatory SNPs.

We identified associations of SNPs spanning over several genes like the transcription factors NFKB1, EGF, and FGF2 spanning 20.3 Mb of 4q24–26 as well as a cluster of genes on 11q13: CCND1, CCS, and GSTP1 associated with a large number of transcripts. These clusters of SNPs were found to be in strong LD (data not shown). When multiple SNPs in the same gene were studied for associations to groups of transcripts, the latter were found often to share a common pathway (Table 9, which is published as supporting information on the PNAS web site). In agreement with Morley et al. (1), who reported 14q32 as a master regulator locus, we found 12 SNPs in that locus containing cancer-related genes like AKT1, MAPK3, CDC42BPB, TNFAIP2, and TRAF3 associated with the expression of 38 genes at P < 0.01.

Genomic Clusters of Associated Transcripts.

Some of the transcripts associated with the same SNP or set of SNPs tended to cluster to common chromosomal regions, as observed in the Centre d’Etude du Polymorphisme Humain families (1). For instance, IL8-associated transcripts KPNA2, FALZ, HN1, and SLC9A3 reside on 17q23-q25; MAPK8-associated transcripts RARRES and SIAT4C reside on 11q23-q23.3, etc. If we plot the P values as a measure of the association between different SNPs and the associated transcripts across all chromosomes we observe a nonrandom distribution, as exemplified in Fig. 8, which is published as supporting information on the PNAS web site, for a stretch of 25 Mb on 17p21. This observation suggests that even in trans interactions are not at random but with transcripts that tend to cluster together on other chromosomes.

Coherent Expression Patterns.

We also studied the pairwise correlations of the expression levels of transcripts associated with the same SNP at P < 0.01 or a set of SNPs for the same gene (Fig. 3). For many such sets of transcripts, the average pairwise correlation values were significantly higher than the random expectation. Correlations were computed as described in Materials and Methods, and the observed values were compared to a distribution obtained by randomly drawing 100 transcript sets of the same size. The observed correlation z scores for associated sets of transcripts, together with properties of the null distributions, are shown in Fig. 3A and Table 3. Furthermore, we performed the same analysis using an independent breast cancer expression data set (Fig. 3B) (17). For several sets of transcripts, the z scores of transcript-to-transcript correlation, as computed from the independent data, were highly significant (Fig. 3 and Table 3), strongly validating the observed associations in our data set. Sets with high z scores in both data sets include the set of transcripts associated with SNPs in ABCB1, BAK, AKT2, and ABCC1 genes. ABCB1 (MDR1) is strongly related to drug resistance in cancer as well as in other conditions. The set of transcripts associated with genetic variants of ABCB1 (Table 10, which is published as supporting information on the PNAS web site) as observed in this study and validated by coherent expression patterns in van’t Veer et al.’s data (17), included an overrepresentation (P < 3E-07) of genes from the Kyoto Encyclopedia of Genes and Genomes human proteasome pathway such as PSMA1, PSMA5, and PSMC6, which might hold the key to better understanding the mechanisms of drug resistance.

Genotype and Haplotype Analyses.

Genotyping was performed as described in Edvardsen et al. (9). LD estimation for the SNPs shown to be associated with transcripts in cis and exemplified in Fig. 5 was performed by using haploview. haploview estimates the maximum-likelihood values of the four gamete frequencies, from which the D′, logarithm of odds, and r² calculations are derived (18). Conformance with Hardy–Weinberg equilibrium is computed by using the Fisher exact test.

Statistical and Bioinformatics Analyses.

Properties of the Association Matrix.

We studied the association matrix P, searching for the following structures:

• Overall overabundance of associated pairs: we assessed the overall significance of the observed association between genotypes and expression phenotypes by comparing it to a null model.

• Potential master regulators: we seek rows and columns in P that have significantly many entries with strong P values. Such structures represent the effect of this locus on many transcripts and suggest the presence of a regulation element. When several SNPs reside in the same gene, these are combined for this analysis. We further study the transcripts that are potentially regulated by a gene or associated with an SNP locus using expression correlation and external information as described below.

• Biclusters: a bicluster is a submatrix of P. We adapt the methods of Tanay et al. (13) to find significantly dense biclusters, in which the participating SNPs share significantly many common associations with the participating transcripts (see Fig. 1B). When several SNP loci, especially when not in LD, associate with the same transcript, this may be evidence of this expression phenotype being a complex trait, affected by several genetic events. When the same set of loci commonly associate with a large number of transcripts, we have cross-confirmation of the individual associations as well as a possible multilocus effect on a pathway or a biological process.

GO and Pathway Analysis.

GO annotations for all transcripts were obtained by using publicly available Biomolecule Naming Service (http://openbns.sourceforge.net), a high-speed directory service that resolves between alias and official gene symbols and links to publicly available databases. We state that a transcript g is linked to a GO term t if t annotates g or is an ancestor of some GO annotation for g. For each GO term t and a list of transcripts L we test the hypothesis that the term t is overrepresented in L against the null hypothesis that distribution of terms is random. This term to transcript–list association is determined by using the hypergeometric distribution [for example, see Benjamini et al. (14)]. Namely, for each term t and list of transcripts L the P value is given by:

where N equals the total number of unique GO annotated transcripts represented in the data set, n equals the total number of unique GO annotated transcripts represented in L, K equals the total number of transcripts linked to the GO term t, and k equals the number of entries in transcript list L linked to the GO term t.

A similar analysis was performed by using pathway information as follows. A pathway database containing 360 curated pathways from various sources such as Kyoto Encyclopedia of Genes and Genomes, BioCarta, and Signaling Pathway Database was used to search for overrepresentation of members of an expanded pathway. Overabundance was assessed by using the hypergeometric distribution similar to the algorithm described above.

Analysis of Correlations of Sets of Associated Transcripts.

The pair-wise correlations between expression of transcripts associated with SNPs in one gene was computed. For each set of SNPs in a gene we considered sets of transcripts associated with an SNP in this gene with P < 0.01. For each set of transcripts, we computed Pearson correlation between each pair of transcripts in the set, and we report the deviation of average correlation in the set from the expected expressed as the z score. We then compared the distribution of these z scores to the expected distribution, computed based on correlations in random sets of transcripts of the same size. The z score was computed as follows: Let C_ij be the correlation of transcripts i and j from the set S. Let C₀ be the average transcript-to-transcript correlation across the entire data set, and let σ₀ be the corresponding standard deviation. For a set of transcripts S with n elements, there are N = n(n − 1)/2 transcript pairs. Then:

Assessment of Statistical Significance.

To assess the statistical significance of our findings, we used the FDR analysis to compare the results to those obtained on random data sets (19). Those were generated by randomly permuting the expression data while leaving the genotype data intact. This randomization process ensures that we keep the structure of dependencies between SNP loci that exist in the original data, as well as between expression vectors. These permuted data were used in assessing overabundance of significant SNP–transcript association pairs as well as in assessing the significance of more complex structures. We used 50 simulations to compute the ANOVA FDR. Genomic intervals on the chromosomes enriched for transcripts with a high average of statistically significant associations for each SNP were estimated as described by Lipson et al. (20).