Patient survival and tumor characteristics associated with CHEK2:p.I157T – findings from the Breast Cancer Association Consortium

Study subjects for survival and pathology analyses

Genotyping

CHEK2:p.I157T was first genotyped by independent studies using various methods including MassARRAY iPLEX Gold (Sequenom, San Diego, CA, USA), TaqMan (Applied Biosystems, Life Technologies, Carlsbad, CA, USA) and Fluidigm (Fluidigm, San Francisco, CA, USA) as listed in Additional file 1: Table S1. Quality control was implemented as follows: each study performed duplicate measurements of at least two samples from each sample plate and genotyped 93 CEPH control DNAs (HAPMAPPT01, Coriell Institute for Medical Research, Cambden, NJ, USA). If a study reported more than two discordant genotyping results of the CEPH DNAs, all genotype data from that study was excluded. Later, p.I157T was genotyped centrally using a custom Illumina iSelect genotyping array for the Collaborative Oncological Gene-environment Study (COGS) [18]. Discordant genotyping results were clarified with Sanger sequencing. CHEK2:c.1100delC was genotyped by independent studies using mainly TaqMan (Additional file 1: Table S1), as described earlier [17].

Pathology analysis

Pathology data was collected from hospital records or from scientific projects within the individual studies, as described previously [19]. Additionally, the TP53 protein expression was measured by individual studies using immunohistochemical staining as described in Additional file 1: Table S3. The pathology data availability and mutation carrier frequencies varied between independent BCAC studies and therefore all analyses were stratified by study. Pathology analyses were performed using R environment for statistical computing version 3.0.2 [20] including packages vcdExtra [21] and meta [22]. Comparisons were made between CHEK2 mutations carriers (heterozygous or homozygous) and non-carriers, for both p.I157T and c.1100delC, as well as between carriers of p.I157T and c.1100delC (Table 1). Associations between the mutations and clinico-pathological characteristics were tested with study-stratified Cochran-Mantel-Haenszel test (mantelhaen.test for categorical characteristics and CMHtest for ordinal characteristics). The category of missing data was not included in these comparisons. Differences in age at diagnosis were tested by meta-analysis of age distribution in independent studies using a random effects model (metacont).

Survival analysis

Survival analyses were performed using the Cox regression [23] as implemented in Stata (Stata/SE 10.1 for Windows, StataCorp LP, College Station, TX, USA) comparing CHEK2 mutation carriers and non-carriers, as described above. Study subjects were considered to become at risk at the time of their first invasive breast cancer diagnosis. The data did not consist entirely of incident cases. Therefore, in order to avoid bias caused by late enrollment, we implemented a method called left censoring, which has been proven to provide robust survival estimates for data, which includes also prevalent cases [24]. Survival analysis endpoints included death of any cause, breast cancer-associated death, distant metastasis relapse, locoregional relapse and second breast cancer. Patients were censored at the end of their follow-up period or at the latest 15 years after the initial breast cancer diagnosis in analyses of overall survival and second breast cancer, but at the latest 10 years in analyses of locoregional or distant relapse-free survival as well as in analyses of breast cancer-specific survival. Patients presenting with distant metastases at diagnosis were excluded from the analyses of locoregional relapse-free survival. All analyses were stratified by study.

In addition to univariate analyses, we performed multivariate analyses, which were stratified by study and age category (≤50 years; >50 and ≤70; >70), and adjusted for tumor grade (1, 2 or 3, ordinal), tumor size (1: maximum diameter less than or equal to 20 mm; 2: more than 20 mm and less than or equal to 50 mm; 3 over 50 mm, ordinal), tumor spread in axillary lymph nodes (0 = negative, 1 = positive) and progesterone receptor (PR) status (0 = negative, 1 = positive). ER was not included in the model, because of the non-linear relationship between tumor ER status and patient survival during the 10 years following the diagnosis; patients with ER-negative tumors have a higher risk of dying from breast cancer during the first 5 years after the diagnosis, but the difference in risk between ER-positive and ER-negative tumors levels out after that period [17, 25]. However, since several studies have reported an association between the two CHEK2 mutations and ER-positive disease [12–14] (Table 1), we performed the survival analyses in a subgroup of patients with ER-positive tumors. Only cases with complete data on the pathological markers were included in the multivariate analyses. Univariate survival analyses were performed also in a subgroup of breast cancer patients with lobular tumors, because of the association between p.I157T and lobular breast cancer [15] (Table 1).

Study subjects for gene expression analysis

Gene expression analyses were performed using a data set of 183 breast tumors from the Helsinki University Hospital (GEO: GSE24450). As described previously, the data set consisted of total RNA samples from 151 tumors from unselected cohorts of breast cancer patients and 32 tumors from additional familial cases hybridized on Illumina HumanHT-12 v3 Expression BeadChips (Illumina Inc., San Diego, CA, USA) [10, 26]. The p.I157T carrier status was defined from peripheral blood samples as described earlier for the BCAC study ‘HEBCS’ (Additional file 1: Table S1). Ten patients were germ line p.I157T carriers and 162 were non-carriers, of which six carried germ line c.1100delC. The c.1100delC carrier tumors were included in the analyses as non-I157T carriers. The p.I157T genotype information was not available for 11 study subjects. These were included in the molecular subtype analysis, but not in differential gene expression or gene set enrichment analysis. The clinico-pathologic characteristics of the 183 tumors are provided in Additional file 1: Table S4.

Gene expression analysis

Gene expression data quality control and quantile normalization was performed in the Bioconductor [27] as described earlier [26]. Data analyses were performed in R version 3.0.2 and Bioconductor packages genefu [28], limma [29, 30] and geneplotter [31]. Probes not mapping to any current Entrez Gene entities (GRCh38.p2) were excluded, resulting in a filtered data set of 20,145 genes.

For determining the intrinsic molecular subtypes, expression data of the fifty PAM50 signature genes was extracted from the filtered data set, median centered and standardized per gene by dividing with the standard deviation of the gene’s expression values. Intrinsic subtypes were defined by Pearson correlation between tumors and the luminal A, luminal B, human epidermal growth factor receptor 2 (Her2)-enriched, basal-like and normal-like centroids as implemented in the genefu package [28, 32]. Hierarchical clustering was performed using the Ward’s method [33]. As a comparison to the subtype classification by gene expression, we used the surrogate clinico-pathologic markers to define the subtypes following the St Gallen 2013 criteria (luminal A: ER+, PR+, Her2-, Ki-67-; luminal B (three marker combinations): ER+, PR-, Her2- or ER+, Her2-, Ki-67+ or ER+, Her2+; basal: ER-, PR-, Her2-; Her2 overexpressing: ER-, PR-, Her2+) [34].

For analysis of differential gene expression the data was filtered by including only genes with highest variation in expression levels over the entire data set (st. dev. ≥ 0.75, 1852 genes). The samples from p.I157T carriers were compared to samples from non-carriers with a moderated t test adjusting for ER, tumor protein 53 (TP53) and Ki-67 protein expression (positive/negative), tumor grade (1, 2, 3, ordinal) as well as histological type (lobular/other). The adjusting covariates were selected from features tabulated in Additional file 1: Table S4 as the most significant factors (p < 0.001) explaining variation in the expression of the 1852 genes as summarized by the first four principal components. Additionally, lobular histologic type was included to avoid bias caused by the association between the p.I157T and lobular type. Data on at least one of the adjusting variables was missing for 12 tumor samples and thus the differential gene expression and gene set enrichment analyses were performed with a set of 160 (ten p.I157T and 150 non-carrier) tumor samples and 1852 genes. Genes with p values below 0.01 were considered to be associated with p.I157T. These were included in a functional enrichment analysis performed using the DAVID functional annotation tool [35]. Functional annotations with Benjamini-Hochberg [36] corrected p values below 0.01 were considered to be significantly enriched.

Gene set enrichment analysis (GSEA) was performed using a java application available at http://software.broadinstitute.org/gsea following the instructions of the user guide [37]. For the GSEA analysis, the 1852 genes were ranked according to a score calculated as the product of log₂(fold change) and log₁₀(p value) from comparisons of p.I157T carrier and non-carrier tumors as described above. All gene sets available at the Molecular Signatures Database (MSigDB) v5.0 [38] were included in the analyses. The p values were corrected for false discovery rate for all other gene sets but the gene sets originating from single publications (‘CGP: chemical and genetic perturbations’ database), which were corrected for the family-wise error rate. Gene sets with corrected p value below 0.05 were considered to be significantly enriched in the p.I157T carrier tumors.