Identification of candidate cancer predisposing variants by performing whole-exome sequencing on index patients from BRCA1 and BRCA2-negative breast cancer families

Background In the majority of familial breast cancer (BC) families, the etiology of the disease remains unresolved. To identify missing BC heritability resulting from relatively rare variants (minor allele frequency ≤ 1%), we have performed whole exome sequencing followed by variant analysis in a virtual panel of 492 cancer-associated genes on BC patients from BRCA1 and BRCA2 negative families with elevated BC risk. Methods BC patients from 54 BRCA1 and BRCA2-negative families with elevated BC risk and 120 matched controls were considered for germline DNA whole exome sequencing. Rare variants identified in the exome and in a virtual panel of cancer-associated genes [492 genes associated with different types of (hereditary) cancer] were compared between BC patients and controls. Nonsense, frame-shift indels and splice-site variants (strong protein-damaging variants, called PDAVs later on) observed in BC patients within the genes of the panel, which we estimated to possess the highest probability to predispose to BC, were further validated using an alternative sequencing procedure. Results Exome- and cancer-associated gene panel-wide variant analysis show that there is no significant difference in the average number of rare variants found in BC patients compared to controls. However, the genes in the cancer-associated gene panel with nonsense variants were more than two-fold over-represented in women with BC and commonly involved in the DNA double-strand break repair process. Approximately 44% (24 of 54) of BC patients harbored 31 PDAVs, of which 11 were novel. These variants were found in genes associated with known or suspected BC predisposition (PALB2, BARD1, CHEK2, RAD51C and FANCA) or in predisposing genes linked to other cancer types but not well-studied in the context of familial BC (EXO1, RECQL4, CCNH, MUS81, TDP1, DCLRE1A, DCLRE1C, PDE11A and RINT1) and genes associated with different hereditary syndromes but not yet clearly associated with familial cancer syndromes (ABCC11, BBS10, CD96, CYP1A1, DHCR7, DNAH11, ESCO2, FLT4, HPS6, MYH8, NME8 and TTC8). Exome-wide, only a few genes appeared to be enriched for PDAVs in the familial BC patients compared to controls. Conclusions We have identified a series of novel candidate BC predisposition variants/genes. These variants/genes should be further investigated in larger cohorts/case-control studies. Other studies including co-segregation analyses in affected families, locus-specific loss of heterozygosity and functional studies should shed further light on their relevance for BC risk. Electronic supplementary material The online version of this article (10.1186/s12885-019-5494-7) contains supplementary material, which is available to authorized users.


Background
Breast cancer is the most common cancer and the leading cause of cancer deaths among women in the world [1]. About 10-20% of all BC patients occur in a familial context, with multiple family members affected across generations [2]. Familial BC susceptibility resulting from deleterious germline variations located on chromosome 17q21 was brought to light through linkage analysis for the first time in 1990 [3]. Since then, many highly penetrant rare variants (with a relative risk of above 10-fold) in BRCA1 (OMIM 113705), BRCA2 (OMIM 600185), TP53 (OMIM 191170), PTEN (OMIM 601728), STK11 (OMIM 602216) and CDH1 (OMIM 192090) to moderately penetrant rare variants (with a relative risk of 2 to 4-fold) in CHEK2 (OMIM 604373), PALB2 (OMIM 610355) [4], BARD1 (OMIM 601593), ATM (OMIM 607585), BRIP1 (OMIM 605882) have been reported. The exact penetrance associated to pathogenic variants in several of these genes is still under investigation. These genes were identified through linkage analysis, positional cloning and/or candidate gene sequencing [5][6][7][8][9][10][11][12][13][14][15][16][17]. Furthermore, with the advent of DNA microarray technology, many low penetrant common variants (with a relative risk often much less than twofold) were unraveled through genome-wide association studies [18]. More recently, thanks to dramatic advances in the speed and scale of nextgeneration sequencing (NGS) technologies combined with sophisticated computation algorithms and a sharp decrease in sequencing cost, a path for the discovery of additional candidate BC predisposing variants has been opened. To name a few, variants in XRCC2 (OMIM 600375), FANCC (OMIM 613899), BLM (OMIM 604610) and PPM1D (OMIM 605100) have been more recently reported as candidate variants with a BC risk through NGS technologies [19][20][21]. The aggregate currently known variants with high, moderate and low penetrance in familial BC susceptibility genes only account for up to 25-50% of all the high-risk BC families. This missing heritability in the remaining 50-75% of BC families [16,22] reflects both the complexity of the BC genetic architecture and the challenges in identifying remaining BC predisposing variants for delivering timely screening, preventive intervention, and precision treatment.
Several studies have revealed that variants in familial BC susceptibility genes like BRCA1, BRCA2, TP53, PALB2, CDH1, PTEN (OMIM 601728), PIK3CA (OMIM 171834), STK11 (OMIM 602216), RINT1 (OMIM 610089) and NF1 (OMIM 613113) are not only associated with BC predisposition, but also with a number of other malignancies [7,9,[23][24][25][26][27][28][29]. In the current study, we hypothesized that in BRCA1 and BRCA2-negative families with elevated BC risk, the analysis of a large array of genes previously associated to (hereditary) cancer syndromes or cancer in general, could likely lead to the identification of additional candidate BC predisposing genes/ variants. Thus, firstly, we identified all rare variants both exome-wide and cancer-associated gene panelwide (492 genes) in 54 BC patients from BRCA1 and BRCA2-negative families with elevated BC risk and compared their relative incidence in 120 geographically matched controls. Secondly, all nonsense, frameshift indels and splice-site variants detected in BC patients within the 492 genes of the panel, which we estimated to possess the highest probability to predispose to BC, were validated on an independent sequencing platform (Roche Junior).

Sample selection
A total of 57 BC patients and 120 controls were considered for this study. Among the BC patients (Additional file 1), 54 were from unrelated BRCA1 and BRCA2-negative families with elevated BC and/or ovarian cancer (OC) risk (i.e. families with two or more affected first-degree relatives) and with a median age at diagnosis of 51 years (range: 36-72). The remaining three BC patients were included as "blinded internal positive controls", each harboring a known germline variant in BRCA1 (NM_007300.3:c.5096G > A), BARD1 (NM_000465.3:c.1921C > T) or PALB2 (NM_024675.3:c. 1571C > G). All geographically matched unrelated controls considered in this study (patients consulted at the same hospital), sequenced according to the same wet lab protocol for cardiac arrhythmias, were unselected for personal or familial history of cancer. The overview of the process of sample preparation, sequencing, analysis and variant validation workflow is presented in ' Additional file 2'. Patient recruitment and blood sampling were performed according to the ethical procedures approved by the institutional ethics committee of the UZ Brussel. Peripheral blood was collected after obtaining a written informed consent for a broad genomic analysis covering also incidental findings in genes predictive for other diseases. Genomic DNA was prepared using Chemagic Magnetic Separation Module I (Chemagen) according to the manufacturer's recommendations.

A virtual panel of cancer-associated genes
After identification of rare variants in whole exomes, we further choose to prioritize variants present in a panel of 492 genes possibly/likely associated with (hereditary) cancer [hereafter called cancer-associated gene panel (CAGP) ( Table 1 and Additional file 3)]. These genes are pooled together from seven gene lists: the well-known cancer susceptibility genes reported by Rahman et al. [30], various BC gene panels reported by Easton et al. [17], genes from BROCA-Cancer Risk Panel (Version 6) [31], Fanconi Anaemia pathway genes reported by Kanchi et al. [32], human DNA repair genes reported by Wood et al. [33], human cancer predisposition genes (GeneRead DNAseq Targeted Panel V2) from Qiagen and genes from the familial cancer database (FaCD, retrieved on 17/02/2015) [34]. Out of the 492 genes from our CAGP, 177 (36%) genes are contributed by at least two gene lists and are mostly known to be cancer susceptibility genes. The remaining 315 (64%) genes are private to a single gene list, mostly from Wood et al. (114 genes) and FaCD (167 genes). Some of these latter genes are not yet clearly associated with (hereditary) cancers (Table 1).

Target-enrichment and next-generation sequencing
For each of the BC patients and controls, one μg of DNA was fragmented using adaptive focused acoustics (Covaris) in order to obtain fragments of approximately 250 base pairs. After DNA end repair and adenylation, oligonucleotides adapters for paired-end sequencing (Illumina) were ligated to both ends of the fragments. Two hundred nanogram of ligated DNA of selected size was PCR amplified and subsequently captured by hybridization for 65 h with the Roche SeqCap EZ Human Exome v3.0 (Roche) Capture Library. After further selection of the targeted fragments through multiple steps of washing, the captured probe-selected DNA was cluster amplified on the Illumina cBot according to manufacturer's protocol (Illumina), using five samples per flow cell lane in order to get sufficient DNA for the subsequent sequencing run. Sequencing was performed on a HiSeq1500 (Illumina) with a paired-end module, generating 125 base reads.

Variant filtration and classification
In-house Python script was used for variant filtration in three steps. Firstly, variants were only retained if they passed VQSLOD (tranche sensitivity threshold of 99.9%) and are located in the exons or at the splice-sites (±2 bp from the exon-intron border). In addition, we required a 10X absolute read depth at the variant position, at least two reads harboring the variant and a variant allele ratio between 20 and 80% along with a minor allele frequency (MAF) ≤1% in any of the population databases (mentioned earlier). Further, we assumed that those variants present in > 10% both in BC patients and controls most likely resulted from sequencing or alignment errors or they should be common variants exclusively in our study population (and thus missed by the MAF restriction). Thus, these variants were removed. Furthermore, missense variants were classified as "probably damaging" (pph2-prob ≥0.957), "possibly damaging" (0.453 ≤ pph2-prob≤0.956), or "benign" (pp2_hdiv ≤0.452) according to PolyPhen-2 (HDIV) [39] in silico prediction scores. Secondly, exome-wide variants that passed all the filters in the first stage were selected for their presence in genes of the CAGP. Lastly, frame-shift indels, nonsense and splice-site variants (hereafter collectively called potentially Protein Damaging Allelic Variants (PDAVs) as they have the highest probability to cause loss of protein function and thus to be associated to BC predisposition) that are present in genes of the CAGP were further validated.  T, TDG, TDP1, TDP2, TELO2, TERC, TGFBR2, TMC6, TMC8, TNF,  TNFRSF11A, TNFRSF13B, TOP3A, TOP3B, TOPBP1, TP63, TREX1, TREX2,  TRIM32, TRPS1, TTC8, TWIST1, TYR, TYRP1,

Variant validation
For validation of the PDAVs obtained from the Illumina platform using capture-based library enrichment system, an orthogonal approach using amplicon-based library enrichment on a 454 platform from Roche (Junior) was performed. Primer pairs were designed in order to amplify DNA fragments (amplicons) that contain the desired variants. One primer of the primer pair was designed towards intronic regions, when possible, to avoid amplification of processed pseudogenes. In addition, BLAST of the target sequence was performed in order to choose only primer pairs that specifically amplify the target region meanwhile avoiding non-specific or pseudo-gene amplification. Furthermore, primers binding to target sequences containing SNPs with a MAF > 1% were avoided. For variant analysis, SeqNext software (JSI medical systems) was used.
Subsequently, we investigated whether an exome-wide enrichment can be observed in the number of variants when comparing BC patients to controls (Student's t-test or Welch's t-test). No significant difference was observed in the average number of variants between BC patients and controls either by pooling all the variant types together (BC patients: controls; 420.81: 426.83, p = 0.3071) or by separately analyzing each sub-type of variants [synonymous SNVs In the next step, we only considered the variants present in the 492 cancer-associated genes from the CAGP panel (see methods). In the BC patients, after filtering, we retained 240 synonymous SNVs, 8 in-frame indels, 13 frame-shift indels, 6 splice-site SNVs, 14 nonsense SNVs and 195 + 94 + 215 missenses SNVs (predicted as "probably damaging", "possibly damaging" and "benign", respectively). In the controls we retained 589 synonymous SNVs, 21 in-frame indels, 23 frame-shift indels, 20 splice-site SNVs, 13 nonsense SNVs and 398 + 174 + 417 missense SNVs (see Additional file 5). When comparing the average number of variants in BC patients versus controls, we observed that the average number of nonsense SNVs was more than twice higher in BC patients [BC patients: controls; 0.26:0.11; ratio = 2.39; p = 0.0287 (0.0688 with Welch correction)], whereas no obvious enrichment could be observed in the other sub-types of variants (see Additional file 6).
To investigate further whether specific genes are more frequently mutated in our BC patients compared to controls, we selected exome wide all the genes harboring high impact mutations (PDAVs) in at least two BC patients (see Additional file 7). Among the 95 genes selected, five (FAM11B, GRAMD2, SP100, USP45 and ZNF534) can be considered candidate BC predisposing genes as they were mutated in three BC patients (out of 54) but not in any of the 120 controls (Additional file 8). Two other good candidate genes are ASPH and C17orf80 as they harbored PDAVs in respectively five and four BC patients and only one control sample (Additional file 8). All PDAVs found in these 7 candidate BC predisposing genes were visually verified using the Integrative Genomics Viewer (IGV) [40].

Validation of PDAVs within the CAGP
PDAVs resulting in dramatic changes in protein structure and function have the highest chance to be associated with BC predisposition. Those PDAVs that were detected in BC patients, passed the filters, and are located within the genes of the CAGP, were further validated on an independent sequencing platform (see methods) and were also reviewed manually using IGV [40]. Thirty-one out of 33 PDAVs present in 24 out of 54 BC patients (~44%) passed the validation step (Table 2), of which 11 PDAVs are not reported in dbSNP147. Among the BC patients with a PDAV, eighteen harbored a PDAV in a single gene, five harbored PDAVs in two genes and one harbored a PDAV in three genes. Furthermore, all five splice-site SNVs (Additional file 9) were considered disruptive by the in silico web-based tool "Human Splicing Finder" [41] (http://www.umd.be/HSF3/, release 3.0). Three of the 26 genes harboring PDAVs in the BC patients were also found mutated in the control samples, suggesting that these genes (ABCC1, BBS10 and PDE11A) are not involved in cancer predisposition (compare Additional files 10 and 11), In addition, the pathogenic variants present in the three internal positive control samples included in this study were also identified.

Discussion
It is expected that exome-wide NGS analysis of a germline DNA sample will reveal many variants when compared to a haploid reference genome, even when only rare variants (MAF ≤ 0.01) are taken into consideration. However, when comparing the total number of variants detected in two individuals of the same ethnicity we do not expect to find significant differences. We confirmed this assumption by (using the same wet bench and dry bench approaches) comparing the average number of variants found in persons belonging to two groups living in the same area (patients recruited in the same hospital): BC patients belonging to elevated risk BC families and controls not selected for personal or familial history of cancer but for cardiac arrhythmias. The ratio of average number of observed (rare) variants in both groups is very close to one for all types of variants ( Fig. 1 (red) and Additional file 6) except for splice site and nonsense variants (0.86 and 0.88, respectively), probably because of the relatively small number of splice site and nonsense variants detected per BC patient and control. When focusing exclusively on the genes of the CAGP, similar observations were obtained (Fig. 1 (blue) and Additional file 6) except for the category of nonsense variants, where more than a two-fold excess of nonsense variants was detected in BC patients compared to controls (ratio = 2.39). Although a larger sample size is a minimal requirement to reach statistical significance, our data suggest that the nonsense variants found in excess in the genes of the CAGP among the BC patients (compared to controls) are implicated in the molecular mechanism modulating BC risk(about 50% of these nonsense variants). If the increased number of nonsense variants seen in BC patients is associated with increased cancer risk, one would expect that these nonsense variants will be more frequently identified in genes functionally correlated with the cancer predisposition process. To verify this assumption, the PANTHER over-representation Test (Released 20,171,205) [42] was used with a false discovery rate (FDR) < 0.05. This over-representation test compares a test gene list to a reference gene list and determines whether a particular class of genes (e.g. those associated to a specific biological process) is overrepresented or underrepresented. We found that genes involved in the DNA repair process, namely inter-strand cross-link repair (FDR: 4.94E-02), double-strand break (DSB) repair via nonhomologous end joining (FDR: 4.35E-02), non-recombinational repair (FDR: 4.42E-02), DSB repair (FDR: 8.35E-02) were overrepresented in BC patients while not in the controls. Indeed, four nonsense variants (out of 14) found in BC patients were found in genes involved in the DSB repair process while only one such variant (out of 13) was found among the controls. It remains unclear for us why the same phenomenon is not observed with the frameshift indels. It is possible that false positive indel calls masked a possible enrichment of the true positive frameshift indels.
Our exome wide analysis revealed only seven genes (ASP, C17orf80, FAM111B, GRAMD2, SP100, USP45 and ZNF534) with high impact mutations (PDAVs) in three (or more) BC patients while comparable mutations were not found (or only once) among the control samples. None of these genes was reported to possess cancer predisposing properties and therefore not included in the CAGP. Gene Ontology (GO) annotation [43] for molecular function, biological processes and Reactome Pathways indicated that SP100 and USP45 are involved in DNA repair while ZNF534 is involved in DNA-templated regulation of transcription, making them good candidate cancer predisposing genes. No molecular function or biological process was annotated to C17orf80 and FAM11B, whereas ASPH and GRAMD2 were reported to be involved in calcium homeostasis and transport (Additional file 8).
When restricting our variant analyses performed on BC patients to the 492 genes of the CAGP, we found novel as well as known PDAVs in several genes known or suspected to be BC predisposing (Table 2 and Additional file 10). Genes participating in DNA DSB repair process e.g. PALB2, BARD1, CHEK2 and RAD51C are particularly intriguing as DSB repair process defective tumors can be selectively targeted by PARP (poly (ADP-ribose) polymerase) inhibitors resulting in synthetic lethality [44][45][46]. We also found PDAVs in genes linked to DNA repair, FA or occurring in some types of cancers but not well studied in the context of familial BC (Table 2 and Additional file 10). These candidate BC predisposing genes are also interesting to scrutinize further in familial BC setting as it is known that familial BC susceptibility genes can also predispose to multiple cancers [30]. Furthermore, we detected PDAVs in genes associated with other hereditary syndromes but not clearly related to cancer ( Table 2 and Additional file 10). These genes are mostly derived from the FaCD panel, which is uncurated. PDAVs detected in the CAGP from control samples but not present in the BC patients are listed in Additional file 11. Only about 44% of the BC patients were found to harbor a PDAV in one (and exceptionally in two or three) gene(s) of the CAGP in this study. Eleven out of 31 PDAVs detected were not reported in dbSNP147 and therefore considered novel. We should keep in mind that these PDAVs are not necessarily BC predisposing. Therefore, their cancer predisposing attributes should be further investigated in much larger cohort /case-control studies or by performing co-segregation analyses in positive families (if sufficient families are available). Although in this study we mainly focused on candidate PDAVs found in genes of the CAGP (which only accounts for 2% of the full exome), we must remain aware that missense variants in genes of the CAGP but also PDAVs and missense variants outside this gene panel may also predispose to BC. For instance, we identified 7 genes not represented in the CAGP in which PDAVs were over-represented in the BC patient cohort. Moreover, BC predisposition may not necessarily rely solely on the presence of one particular variant in the family but may result from combinatorial interactions between several variants. Indeed, it has been proposed that in the majority of BC families, BC predisposition could be polygenic in nature and the contribution of several variants located in genes associated to moderate or low risk could be responsible for the increased susceptibility to BC [47]. The mechanism how these different variants cooperate at the molecular level to create an increased BC risk is a matter of further investigation [48].

Conclusions
On average, twice more nonsense variants were found in BC patients than in controls when analyzing the genes from the CAGP. Moreover, GO analysis (biological process) of the genes accumulating those nonsense variants indicated that genes involved in the DSB repair process were overrepresented in the BC patients but not in controls. Comparable observations were not made for the other variant types in the CAGP, nor when considering the whole exome. Taken together, our observations might indicate that a nonsense variant found in the CAGP of a BC patient has more than 50% chance to be associated with BC risk while similar conclusions cannot be drawn for "probably/possibly damaging" missense or frameshift mutations. Larger case-control studies should be performed to confirm these assumptions and validate our candidates. This preliminary study in 54 BC patients from BRCA1 and BRCA2-negative BC families with elevated cancer risk identified candidate BC predisposing PDAVs (known as well as unknown) in 30 genes; PALB2, BARD1, CHEK2, RAD51C, FANCA, RINT1, EXO1, RECQL4, CCNH, MUS81, TDP1,