Expression Quantitative Trait loci (QTL) in tumor adjacent normal breast tissue and breast tumor tissue

Adequacy of FFPE expression quantification

It is understood that expression analysis of FFPE samples can be suboptimal owing to RNA degradation and other chemical and structural effects of fixation [10]. To build confidence in our expression quantifications, we acquired RNA-seq expression data on 1020 fresh-frozen breast cancer samples as archived at TCGA with date tag 2015-11-01. HUGO gene identifiers were matched between NHS and TCGA platforms for 9316 expression quantification targets. We focused on the capacity of expression measures to discriminate tumor hormone receptor status, as all samples were classified as estrogen receptor positive or negative (ER+ or ER-), and progesterone receptor positive or negative (PR+ or PR-). The 2 x 9316 F-statistics for the 3 d.f. tests of equal mean expression over ER+, ER-, PR+, PR- samples exhibited positive rank correlation (Spearman r = 0.43, p < 10^-10). A more focused assessment, based on 17 elements of the 21-gene Oncotype DX recurrence signature common to the NHS and TCGA platforms, yielded a Spearman's r = 0.85, p < 10^-10; see Fig 2. We conclude from this analysis that the RNA extracted from the NHS FFPE samples is suitable for differential expression analysis.

Overall eQTL results

Seventy-one (71) breast cancer risk SNP were analyzed in this report; see S1 Table for their identifiers and information on genomic context. Separate analyses were performed for ER+, ER-, tumor adjacent normal, and contrast between tumor and tumor adjacent normal tissue. Test results were filtered to a maximum false discovery rate (FDR) of 10% per SNP, per tissue type. With this procedure we identified, using an additive genetic model based on called genotypes or imputed allelic dosage, a total of 27 trans-acting SNPs, exhibiting effects on expression of 369 unique Affymetrix Human Transcriptome Array (HTA) version 1.0 transcript clusters corresponding to 217 annotated genes. All but 10 significant associations involved SNPs and genes on different chromosomes; for the significant same-chromosome associations, the minimum SNP-gene separation was 6.4Mb. Consequently all associations are referred to as "trans-eQTL". Table 1 enumerates SNPs with significant trans-eQTL in terms of location, tissue-type analyzed, structural context of variants within genes, or epigenetic contexts of intergenic variants (as estimated with ChromHMM (11) on the HMEC cell line), the number of trans associations identified in different tissue types, minor allele frequency, and, when available, CADD score for variant deleteriousness, reported on the PHRED scale[11]. Findings in Table 1 are presented in the order: tumor vs. normal directly compared (5 SNP, column TUM/NOR), ER+ (4 SNP), ER- (7 SNP), multiple tissue findings (9 SNP), findings in normal tissue only (2 SNP). Table 1 condenses information to one record per SNP; full information on all significant associations is provided in S2 Table.

We have created a browser-driven utility (bceBrowse, in Bioconductor package bceQTL, forthcoming, introductory video at https://www.youtube.com/watch?v=fJDXI5M7_mQ) that allows searching, sorting and plotting the FDR-filtered association results to facilitate tabulating and navigating of our data. A Circos display of trans-eQTL relationships for expression measured in ER- tumors is provided in Fig 3. Supplementary S1 and S2 Figs depict the relationships for findings in paired tumor/adjacent normal, and ER+ tumors.

To illustrate variation in expression present in trans associations inferred in this study, the bceBrowse utility will display boxplots for selected SNP-gene pairs. Fig 4 provides examples of 4 trans-eQTL findings, each stratified by ER status and tissue type (tumor vs. adjacent normal).

Functional QTLs

We tested SNPs for association with gene groups defined by the Molecular Function (MF) Gene Ontology terms to identify pathway-level associations for the 71 breast cancer loci. The associations presented in Table 2 demonstrate that a unit increase in the dosage of a SNP modifies the mean expression of the genes in a pathway-specific gene set categorized as the MFs in GO.

Study population

For any report of breast cancer, written permission was obtained from participants to review their medical records to confirm the diagnosis and to classify cancers as in situ or invasive, by histological type, size, stage of disease and presence or absence of metastases. The Human Subjects Committee at Partners Healthcare System and Brigham and Women’s Hospital in Boston, Massachusetts, have reviewed and approved this study.

The Nurses’ Health Study cohort was established in 1976 when 121,701 female US registered nurses ages 30–55 responded to a mailed questionnaire inquiring about risk factors for breast cancer. Study participants reported new breast cancer diagnoses on subsequent biennial questionnaires. Pathology reports were also reviewed to obtain information on estrogen receptor (ER), progesterone receptor (PR), and Human Epidermal growth factor Receptor 2 (HER2) status. Collection procedures for the breast cancer tissue blocks have been described in detail previously[34]. Participants complete a biennial questionnaire on dietary and lifestyle factors, reproductive factors, anthropomorphic measures, medication use, and health outcomes. Health outcomes include nonfatal incident diseases such as cancer. The follow-up rate in this cohort has been over 90%[35].

The Cancer Genetic Markers of Susceptibility (CGEMS) project (http://www.cgems.cancer.gov) is an NCI initiative to conduct genome-wide association studies (GWAS) to identify genes involved in breast and prostate cancer[36]. As a part of this project 1,145 cases were genotyped using the Illumina 550 array (data is publically available in dbGAP). For the current analysis, we identified invasive postmenopausal breast cancer cases diagnosed from 1990–2004 with GWAS data from the CGEMS project and sufficient RNA for expression profiling in breast tumor and tumor adjacent normal breast tissue.

Informed consent was obtained from all Nurses’ Health Study participants, including consent to collect and use blood samples and tissue specimens for genetic and genomic research[1,35]. This study was reviewed by the Internal Review Board and approved by the Human Subjects Committee at Partners Healthcare System and Brigham and Women's Hospital in Boston, Massachusetts.

SNP selection

The 71 germline genetic variants used in the eQTL and fQTL analyses were the published loci from breast cancer GWAS and the U19/GAME-ON meta-analysis [1–3] (http://gameon.dfci.harvard.edu). Genotyped and imputed SNP dosages were obtained from the CGEMS study[1] breast cancer cases.

Gene expression assay methods

RNA was extracted from multiple 1 or 1.5 mm cores (median number of cores was 5 in normal samples and 3 in tumor samples) taken from representative portions of tumor or adjacent normal from the FFPE blocks. The RNA extraction was performed using the Qiagen AllPrep RNA isolation kit. Total RNA was amplified using the NuGEN WT-Ovation FFPE System (NuGEN, San Carlos, CA). Total RNA was used to synthesize double-stranded complementary DNA (cDNA) utilizing a random priming method and then double-stranded cDNA was fragmented, labeled and hybridized to the Affymetrix Glue Grant Human Transcriptome Array (HTA 1.0) pre-release version from Affymetrix, Santa Clara, CA[37]. CDNA yield was measured in subset of samples prior to hybridization as a quality control measure. After hybridization and washing, processed slides were scanned with the Affymetrix GeneChip Scanner 3000 7G at the Dana-Farber Cancer Institute Microarrray Core Facility.

The Glue Grant Human Transcriptome Array was developed for high-throughput clinical studies, allowing for comprehensive examination of gene expression and genome-wide identification of alternative splicing as well as detection of coding SNPs and noncoding transcripts(23). The microarray can detect 34,834 lncRNAs and 39,224 transcript clusters (genes) curated from the most authoritative databases such as RefSeq and Ensembl as well as the literature. Repeat sequences and ncRNAs shorter than 200 bp were deleted. Each transcript was represented by 119 unique probes (on average) to improve statistical confidence. Each transcript was represented by a specific exon or splice junction probe that can accurately identify individual transcripts.

Quality control analysis

The quality control procedures include patient sample filtering, probe filtering, PCA to adjust for batch-to-batch variation in sample preparation (e.g. RNA extraction), and evaluation of subtype concordance by immunohistochemistry and HTA probe. In Fig 1, we summarize the number of NHS participants at each stage of the study, from 774 eligible breast cancer specimens to the final sample size (639 specimens with gene expression and genotype data) included in the analysis. To ensure data quality, we conducted quality control analyses (correlation of four independent technical replicates across nine assay plates r≥0.93), biological data checks (probe concordance with markers measured by immunohistochemstry), as well as probe filtering and Principal Components Analysis (PCA) to correct for assay batch-to-batch variation. We included four independent breast tumor samples as technical replicates in each assay plate. The technical replicate breast tumor specimens were identified from Beth Israel Deaconess Medical Center. Sufficient RNA was extracted using the Qiagen AllPrep Kit at a single time point to include in all assay plates. First, we excluded ten specimens (one cel file per specimen) that could not be read by hGlueQC. We normalized the data using the robust median average (RMA) method. Next, we evaluated sample quality and excluded 124 NHS specimens with less than 0.55 area under the curve (AUC) threshold for distinguishing a positive and negative probe signal. Then, we evaluated the association of the first fifty principal components (PCs) in the combined tumor adjacent normal and tumor specimens after robust median average normalization with assay batch. We clipped the first PC associated assay batch

Adequacy of FFPE-based expression assays

The NHS provides expression data on 326 FFPE samples. To assess the biological interpretability of the NHS expression assays based on FFPE tissue, we obtained the RSEM-based RNA-seq gene quantifications on 1020 non-FFPE samples from the TCGA BRCA cohort. Using HUGO identifiers for genes, the two platforms had 9316 features in common. For platforms p∈[22] we fit (using regularized linear models [38]) the 9316 models

with g indexing genes and i indexing individuals assayed on platform p. Here y_igp denotes the normalized expression level for gene g measured on subject i in platform p, x_i is a dummy 3-vector indicating whether patient i is +-, -+, or ++ for ER/PR tumor status, and e_igp is a homoskedastic error assumed independent across individuals. The (regularized) F statistics for H_o:β_gp = 0 were computed for each platform as a basic measure of assay sensitivity to differences in tumor biology. Owing to the different sample sizes, the denominator degrees of freedom for the F statistics are very different (322 for NHS and 1016 for TCGA). We use the statistics to rank the models, and thereby genes, according to their capacity to distinguish expression variation between samples derived from tumors of different hormone receptor types.

CADD annotation

Combined Annotation Dependent Depletion (CADD)[11] is a tool for scoring the deleteriousness of SNPs (included in Table 1). CADD integrates multiple annotations into one metric by contrasting variants that survived natural selection with simulated mutations. C-scores are highly associated with allelic diversity, pathogenicity of both coding and non-coding variants and experimentally measured regulatory effects and also highly rank causal variants within individual genome sequences. CADD can quantitatively prioritize functional, deleterious and disease causal variants across a wide range of functional categories, effect sizes, and genetic architectures.

Statistical analysis

To identify regulatory variants, we investigated the association of 71 breast cancer risk alleles (each entered as a single continuous dosage variable as counts of the minor allele) in linear regression models to identify associations with normalized gene expression levels (eQTL) or biological pathways (fQTL). Covariate information for the breast cancer cases was obtained from biennial questionnaires or extracted from medical records or a supplemental questionnaire. We included patient’s age at diagnosis and year of diagnosis along with array plate identifier as covariates in the multivariate linear regression model, and assumed an additive effect of the SNP on the gene expression measurement. For s = 1,…, 71 enumerating the risk SNP, g = 1,…, 26001 HTA 1.0 transcript clusters, and individuals i contributing tissues of a given type, let y_ig denote the expression level of gene g for subject i, $X_i^t$ denote the vector of subject's age of diagnosis, year of diagnosis, array plate number, and contribution to the principal component adjustment for expression heterogeneity, and let D_is denote the dosage of the risk allele for SNP s. Statistical significance of the association between SNP s and gene g is determined using the FDR based on the moderated t-statistic for testing the null hypothesis Ho: β_2sg = 0 in the models

where ∈_isg denotes a random quantity with mean zero and finite variance ${\sigma }_{sg}^2$. For the analysis of paired tumor-adjacent normal samples, the model is elaborated to include fixed effects of tissue type (dichotomous, tumor or normal) and allelic dose, and the parameter of interest is the coefficient of the product of tissue type code and allelic dose.

Association tests were conducted separately for breast tumor (total tumor n = 376; by estrogen receptor 1 alpha [ERα] by IHC: ER+ tumor n = 262, ER- tumor n = 70) and tumor adjacent normal samples (n = 264) using Affymetrix Power Tools and the Bioconductor package linear models for microarray data (LIMMA)[38]. LIMMA computes modified gene-specific t-statistics for general linear predictors of mean expression, and incorporates the Benjamin-Hochberg FDR multiple testing correction technique. Among 250 paired tumor-normal specimens, we also conducted regression analysis testing for an interaction between SNP genotype and tissue type, adjusting for all covariates noted above and including a fixed effect of individual contributing the paired samples.

fQTL analysis

Functional QTL (fQTL) analysis tests the hypothesis that SNPs influence not only the expression of single genes, but multi-gene processes such as those represented in pathways and annotated in pathway databases such as Reactome[39,40] or the Gene Ontology: classification system. Because we are testing for association between SNPs and groups of genes, fQTL analysis identifies significant trans-associations between sequence variants and biological mechanisms or pathways and treats the measurements of gene expression measures on individual genes as repeated measures on processes. This procedure could be regarded as treating the expression data for all genes in a group as a multivariate response. This solution is related to but not identical to approaches for gene set testing in transcriptomics. The novel question we are attempting to answer is whether there are particular genetic variants that are associated with the expression of a functionally-related group of genes.

The rationale underlying the fQTL model is that a locus or loci within the genome may affect the activity of an intermediary genomic factor, such as a transcription factor, micro RNA, or gene promoter region and that this intermediate element may then affect the expression of a set of genes.[14,22–29] Procedures confronting this generic concern include methods by Wolfringer[41], Goeman's global test [42], Barry's SAFE[43], GSEA[44] and several procedures due to Wu and colleagues [45,46], We defined “pathways” using the GO Molecular Function classifications and then were standardized using the mean and standard deviation values for those genes, which were used to calculate Z-scores. Each pathway is then represented as a vector and the size of this vector depends on the number of samples in the analysis. We construct and fit a linear model of the form:

where YrDx is the year of diagnosis, AgeDx is age of diagnosis and i is the gene set, j is for the gene within the gene set and k is for the sample. The risk parameter β₃ represents the increase in the mean pathway specific Z-score for each unit increase in the SNP-dosage. It is important to note that this model assumes that genes within a biological class have the same variance, and that, by treating genes within a class as repeated measures on a process, we increase the power of the analysis by reducing the effective number of tests in comparison to a more standard eQTL analysis.

In our fQTL analysis we tested for association between SNPs and the expression of 396 gene sets defined by the Gene Ontology Molecular Function classification for both the 376 tumor and the 263 tumor adjacent normal samples.

Multiple comparisons

To adjust for multiple comparisons, we used the False Discovery Rate (FDR) for both eQTL and fQTL analyses and present results for FDR < 10%[47]. The eQTL analyses incorporated the Benjamini-Hochberg FDR multiple hypotheses correction technique[38]. For the fQTL analyses, we calculated the FDR through a permutation procedure based on the genotype information[48]. By permuting only the genotype information of the samples, the correlation structure within gene expression and genotype data is preserved.