Exome sequencing identifies novel susceptibility genes and defines the contribution of coding variants to breast cancer risk

UK biobank

The UK Biobank is a population-based prospective cohort study of more than 500,000 samples. More detailed information on the UK Biobank is given elsewhere^17,18. WES data for 200,000 samples were released in October 2020¹⁹. QC metrics were applied to Variant Call Format files¹⁹. At the genotype level, SNPs were excluded with sequencing depth <7 or heterozygous allele balance <0.15 or >0.85. Indels were excluded with sequencing depth <10 or allele balance <0.2 or >0.8. On males’ X chromosomes, depth filters were reduced to 5 for SNPs and 7 for indels. Samples with missing calls for >15% and variants with missing calls for >15% of samples were excluded. Variants with Hardy Weinberg Equilibrium p-value 10^-15 were also removed.

Samples with disagreement between genetically determined and self-reported sex, sex aneuploidy, or excess relatives in the dataset were excluded. Excess relatives were identified by considering pairs of individuals with kinship>0.17. If one individual in a pair was a case and one was a control then the case was preferentially selected; otherwise, one individual was selected at random. Genetic ancestry was estimated using genetic principal components and the Gilbert-Johnson-Keerthi distance algorithm²⁰. If genetic principal components were not available, self-reported ethnicity was used. Samples of ancestry other than European were excluded. The final dataset for analysis included 181,992 samples with 100,068 Females. Cases were defined by having invasive breast cancer (ICD-10 C50) or carcinoma-in-situ (D05), as determined by linkage to national cancer registration (NCRAS), or self-reported breast cancer. Both prevalent and incident cases were included. Only breast cancers which were an individual’s first or second diagnosed cancer were included as cases. By this definition, 8,043 female and 45 male cases were included.

The Breast Cancer Association Consortium datasets

The BRIDGES and PERSPECTIVE samples were from studies in the BCAC (BRIDGES: 8 studies, PERSPECTIVE: 3 studies; Supplementary Table 1). Most samples were previously included in a targeted panel sequencing project¹. Phenotype data were based on the BCAC database v14. Samples were oversampled for early-onset (age of diagnosis below 50 years) or family history of breast cancer. Cases were preferentially selected to have information on tumor pathology. Samples with previously identified pathogenic mutations in BRCA1, BRCA2 or PALB2 were excluded.

For BRIDGES, library preparation was conducted in the 3 laboratories using the Nextera DNA Exome kit (Illumina) for tagmentation, barcoding and amplification steps. Subsequently, 500ng of DNA per sample were pooled in 12-plex and concentrated using a vacuum system. Afterwards, hybridization capture reagents for DNA libraries were used for overnight hybridization with the xGen Exome Research panel (Integrated DNA technologies, IDT), capture and amplification. Barcoded pooled libraries of 96 samples were sequenced on each lane of a NovaSeq 6000 S4 flowcell (Illumina) using NovaSeq Xp 4-Lane Kit (2×100bp).

For PERSPECTIVE, library preparation was conducted using Agilent SureSelect Human all exon V7 (48.2Mb). Barcoded libraries of 88 samples were sequenced on a NovaSeq 6000 S4 flowcell (Illumina) using NovaSeq Xp4-lane Kit (2×100bp).

The same pipeline for variant calling was applied to both the BRIDGES and PERSPECTIVE data and followed the GATK (Genome Analysis Toolkit) best practices. Briefly, raw sequence data (FASTQ format) were pre-processed to produce BAM files. This involved alignment to the reference genome, identification and removal of duplicate read pairs from the same DNA fragments, and base recalibration. The base recalibration included the generation of a base quality score recalibration table, later applied to the read bases to adjust their quality scores and increase the accuracy of the variant calling algorithms. An intermediate and informal Quality Control was performed for a sanity check, including coverage and alignment mapping metrics. Variants were then called using Haplotype Caller for the whole exome. This was later split into chromosomes for the analysis due to file size constraints. Variants were further filtered using Variant Quality Score Recalibration (VQSR), a proprietary algorithm from GATK that applied new calibration scores independently at SNP and indel variant levels.

From the final dataset, samples and variants were excluded based on coverage, allelic balance, and Hardy-Weinberg equilibrium, using the same filtering as for UK Biobank. We also excluded samples where the genotypes were inconsistent with previous array genotyping or targeted sequencing data^1,2.

The BRIDGES study sequenced 6,912 samples, of which 3,461 cases and 3,200 controls remained in the final dataset after QC. The PERSPECTIVE study sequenced 10,523 samples, of which 4,777 cases and 5,210 remained in the final dataset.

Data preparation

For both the UK Biobank and BCAC datasets, Ensembl Variant Effect Predictor (VEP) was used to annotate variants²¹. Annotations included the 1000 genomes phase 3 allele frequency, sequence ontology variant consequences and exon/intron number. For each gene, the MANE Select²² transcript was used if it was available for that gene, or the RefSeq Select transcript²³. Annotation files were used to identify PTVs and rare (allele frequency <0.001 in both the 1000 genomes dataset and the current dataset) missense variants. PTVs in the last exon of each gene were excluded as these are generally predicted to escape Nonsense-Mediated mRNA Decay (NMD). VEP was also used to annotate missense variants by Combined Annotation Dependent Depletion score (CADD)²⁴. Here CADD ≥20 was used to define variants predicted to be deleterious.

Burden Test analysis

Association analyses were carried out for each gene separately for PTVs, rare missense variants, and predicted deleterious rare missense variants (defined by CADD score ≥20) and PTVs combined. The main association analyses were burden tests in which genotypes were collapsed to a 0/1 variable based on whether samples carried a variant of the given class. That is, and 0 if where g_ij = 0, 1, 2 is the number of minor alleles observed for sample i at variant j, and p is the number of variants in the gene (thus, heterozygous and homozygous carriers were combined).

Different statistical tests were considered by simulation, and the best method was chosen for each dataset. For the BCAC datasets, we used an exact conditional test^25,26, stratified by country and library preparation method (BRIDGES vs Perspective). Ethnicity was not adjusted for as within each country ethnicity was constant. This method had greater power than the Mantel-Haenszel test²⁷ and Wald test from logistic regression²⁸, and the type-1 error rate was closest to the specified significance level. The exact conditional test was also used for case-control analyses for subtypes and case-only analysis when comparing subtypes e.g., oestrogen receptor status.

Two genes with very common PTVs, NUDT11 and ZNF598, were excluded because missing genotypes (which were treated as non-carriers) led to spurious associations even though the variants passed QC filtering. AFF1 was also removed as the PTV frequency was high in PERSPECTIVE but rare in BRIDGES and UKB. This was likely due to a single PTV artefact within the PERSPECTIVE dataset.

For UK Biobank, we used logistic regression analysis but also incorporated family history as a surrogate for disease status. This markedly improves power since susceptibility variants will also be associated with family history; in particular, it allows information on males in the cohort with a family history of female breast cancer to be utilised. To do this, we treated genotype (0/1) as the dependent variable and family history weighted disease status as the covariate; the latter is defined as d+1/2f, where d=0,1 was the disease status of the genotyped individual and f=0 or 1 according to whether the individual reported a positive first-degree family history. The rationale for this weighting is that, for small effect sizes, the log(odds ratio) associated with a positive first-degree relative is approximately ½ that associated with the disease. All analyses adjusted for the first 10 principal components. For genes on chromosome X, only females were used in the analysis.

To combine the results from the BCAC and UK Biobank datasets in a meta-analysis, the association tests for each gene were converted to Z-scores. The combined Z-score was defined as: . Here, Z_M is the combined z-score, Z_j is the z-score for study j and w_j is the weight associated with study j.

A standard meta-analysis would define the weights w_j using inverse variance or effective sample sizes. However, the effect sizes from the BCAC and UK Biobank may not be comparable, since the BCAC studies oversampled for family history and early age at onset, which may have increased the estimated effect. Furthermore, the UK biobank analysis incorporated family history, which changes the effective sample size as well as the effect estimate. Therefore, we defined weights by using the associations in the known risk gene CHEK2 as a standard: we rationalised that the CHEK2 PTVs provided the best standard as the association is well established and the odds ratio is highly reproducible^1,29,30. Moreover, the odds ratio (∼2-2.5) was representative of the size of effects we hoped to detect for other genes. Thus, we defined , where is the ratio of z-scores for CHEK2 for the BCAC dataset and UK Biobank. A similar approach was applied for the meta-analysis of the other variant categories. Z_M-scores were plotted in Manhattan plots and associated P-values plotted in quantile-quantile plots. The lambda statistic for inflation in the test statistics (based on the median chi-squared statistic) was 0.638 for UK Biobank, 0.549 for BCAC and 0.571 for the meta-analysis, indicating that the tests were somewhat conservative on average.

To investigate the joint effect of PTVs in MAP3K1 and common susceptibility variants in the region identified through GWAS, we accessed imputed genotype data from UK Biobank for the lead SNPs as identified through previous fine-mapping analyses^10,11 We fitted logistic regression models including covariates for PTVs and the lead SNPs and compared the fit of the model, and effect sizes, with the model in which the PTVs or the lead SNPs were excluded.

Data on clinicopathological characteristics of cases in the BCAC dataset was also accessed and the proportion of individuals with specific pathologic features e.g., stage and grade were compared between carriers of variants in a specific gene, e.g., MAP3K1 PTV carriers, and the overall dataset.

Contribution of PTVs to the FRR

We estimated the overall contribution of PTVs to the familial relative risk (FRR) of breast cancer using an empirical Bayesian approach. Given the aggregate frequency p_j of PTVs in a gene is rare, and all PTVs confer the same relative risk , the FRR due to one gene, given p_j and , is: Under the additional assumption that the risks conferred by variants in different genes are additive, the total contribution over J genes is given by: We assumed a prior distribution for effect sizes (log-odds ratios) in which a proportion α of genes are risk associated, and the estimated log-odds ratios β for associated genes have an exponential distribution with parameter η; this distribution was chosen since the distribution of effect sizes is likely to be skewed, with only a small number of genes have a large effect size and most undiscovered genes having smaller effect sizes. An approximate likelihood of the observed carrier count data, by gene, was derived, summed over all genes and maximised numerically to estimate α and η, and hence posterior effect size distributions given the data. The total contribution to the FRR was estimated by integrating the FRR estimates given β_j over the posterior distribution. Further details are given in Supplementary Methods.