A unified analytic framework for prioritization of non-coding variants of uncertain significance in heritable breast and ovarian cancer

Design of tiled capture array for HBOC gene panel

Nucleic acid hybridization capture reagents designed from genomic sequences generally avoid repetitive sequence content to avoid cross hybridization [77]. Complete gene sequences harbor numerous repetitive sequences, and an excess of denatured C₀t-1 DNA is usually added to hybridization to prevent inclusion of these sequences [78]. RepeatMasker software completely masks all repetitive and low-complexity sequences [79]. We increased sequence coverage in complete genes with capture probes by enriching for both single-copy and divergent repeat (>30 % divergence) regions, such that, under the correct hybridization and wash conditions, all probes hybridize only to their correct genomic locations [77]. This step was incorporated into a modified version of Gnirke and colleagues’ (2009) in-solution hybridization enrichment protocol, in which the majority of library preparation, pull-down, and wash steps were automated using a BioMek® FXP Automation Workstation (Beckman Coulter, Mississauga, Canada) [80].

Genes ATM (RefSeq: NM_000051.3, NP_000042.3), BRCA1 (RefSeq: NM_007294.3, NP_009225.1), BRCA2 (RefSeq: NM_000059.3, NP_000050.2), CDH1 (RefSeq: NM_004360.3, NP_004351.1), CHEK2 (RefSeq: NM_145862.2, NP_665861.1), PALB2 (RefSeq: NM_024675.3, NP_078951.2), and TP53 (RefSeq: NM_000546.5, NP_000537.3) were selected for capture probe design by targeting single copy or highly divergent repeat regions (spanning 10 kb up- and downstream of each gene relative to the most upstream first exon and most downstream final exon in RefSeq) using an ab initio approach [77]. If a region was excluded by ab initio but lacked a conserved repeat element (i.e. divergence > 30 %) [79], the region was added back into the probe-design sequence file. Probe sequences were selected using PICKY 2.2 software [81]. These probes were used in solution hybridization to capture our target sequences, followed by NGS on an Illumina Genome Analyzer IIx (Additional file 1: Methods).

Genomic sequences from both strands were captured using overlapping oligonucleotide sequence designs covering 342,075 nt among the 7 genes (Fig. 1). In total, 11,841 oligonucleotides were synthesized from the transcribed strand consisting of the complete, single copy coding, and flanking regions of ATM (3513), BRCA1 (1587), BRCA2 (2386), CDH1 (1867), CHEK2 (889), PALB2 (811), and TP53 (788). Additionally, 11,828 antisense strand oligos were synthesized (3497 ATM, 1591 BRCA1, 2395 BRCA2, 1860 CDH1, 883 CHEK2, 826 PALB2, and 776 TP53). Any intronic or intergenic regions without probe coverage are most likely due to the presence of conserved repetitive elements or other paralogous sequences.

For regions lacking probe coverage (of ≥ 10 nt, N = 141; 8 in ATM, 26 in BRCA1, 10 in BRCA2, 29 in CDH1, 36 in CHEK2, 15 in PALB2, and 17 in TP53), probes were selected based on predicted T_ms similar to other probes, limited alignment to other sequences in the transcriptome (<10 times), and avoidance of stable, base-paired 2° structures (with unaFOLD) [82, 83]. The average coverage of these sequenced regions was 14.1–24.9 % lower than other probe sets, indicating that capture was less efficient, though still successful.

HBOC samples for oligo capture and high-throughput sequencing

Genomic DNA from 102 patients previously tested for inherited breast/ovarian cancer without evidence of a predisposing genetic mutation, was obtained from the Molecular Genetics Laboratory (MGL) at the London Health Sciences Centre in London, Ontario, Canada. Patients qualified for genetic susceptibility testing as determined by the Ontario Ministry of Health and Long-Term Care BRCA1 and BRCA2 genetic testing criteria [84] (see Additional file 2). The University of Western Ontario research ethics board (REB) approved this anonymized study of these individuals to evaluate the analytical methods presented here. BRCA1 and BRCA2 were previously analyzed by Protein Truncation Test (PTT) and Multiplex Ligation-dependent Probe Amplification (MLPA). The exons of several patients (N = 14) had also been Sanger sequenced. No pathogenic sequence change was found in any of these individuals. In addition, one patient with a known pathogenic BRCA variant was re-sequenced by NGS as a positive control.

Sequence alignment and variant calling

Variant analysis involved the steps of detection, filtering, IT-based and coding sequence analysis, and prioritization (Fig. 2). Sequencing data were demultiplexed and aligned to the specific chromosomes of our sequenced genes (hg19) using both CASAVA (Consensus Assessment of Sequencing and Variation; v1.8.2) [85] and CRAC (Complex Reads Analysis and Classification; v1.3.0) [86] software. Alignments were prepared for variant calling using Picard [87] and variant calling was performed on both versions of the aligned sequences using the UnifiedGenotyper tool in the Genome Analysis Toolkit (GATK) [88]. We used the recommended minimum phred base quality score of 30, and results were exported in variant call format (VCF; v4.1). A software program was developed to exclude variants called outside of targeted capture regions and those with quality scores < 50. Variants flagged by bioinformatic analysis (described below) were also assessed by manually inspecting the reads in the region using the Integrative Genomics Viewer (IGV; version 2.3) [89, 90] to note and eliminate obvious false positives (i.e. variant called due to polyhomonucleotide run dephasing, or PCR duplicates that were not eliminated by Picard). Finally, common variants (≥1 % allele frequency based on dbSNP 142 or > 5 individuals in our study cohort) were not prioritized.

IT-based variant analysis

All variants were analyzed using the Shannon Human Splicing Mutation Pipeline, a genome-scale variant analysis program that predicts the effects of variants on mRNA splicing [91, 92]. Variants were flagged based on criteria reported in Shirley et al. (2013): weakened natural site ≥ 1.0 bits, or strengthened cryptic site (within 300 nt of the nearest exon) where cryptic site strength is equivalent or greater than the nearest natural site of the same phase [91]. The effects of flagged variants were further analyzed in detail using the Automated Splice Site and Exon Definition Analysis (ASSEDA) server [38].

Exonic variants and those found within 500 nt of an exon were assessed for their effects, if any, on SRFBSs [38]. Sequence logos for splicing regulatory factors (SRFs) (SRSF1, SRSF2, SRSF5, SRSF6, hnRNPH, hnRNPA1, ELAVL1, TIA1, and PTB) and their R _sequence values (the mean information content [93]) are provided in Caminsky et al. (2015) [36]. Because these motifs occur frequently in unspliced transcripts, only variants with large information changes were flagged, notably those with a) ≥ 4.0 bit decrease, i.e. at least a 16-fold reduction in binding site affinity, with R _i,initial  ≥ R _sequence for the particular factor analyzed, or b) ≥ 4.0 bit increase in a site where R _i,final  ≥ 0 bits. ASSEDA was used to calculate R _i,total , with the option selected to include the given SRF in the calculation. Variants decreasing R _i,total by < 3.0 bits (i.e. 8-fold) were predicted to potentially have benign effects on expression, and were not considered further.

Activation of pseudoexons through creating/strengthening of an intronic cryptic SS was also assessed [94]. Changes in intronic cryptic sites, where ΔR _i  > 1 bit and R _i,final  ≥ (R _sequence – 1 standard deviation [S.D.] of R _sequence ), were identified. A pseudoexon was predicted if a pre-existing cryptic site of opposite polarity (with R _i  > [R _sequence - 1 S.D.]) and in the proper orientation for formation of exons between 10–250 nt in length was present. In addition, the minimum intronic distance between the pseudoexon and either adjacent natural exon was 100 nt. The acceptor site of the pseudoexon was also required to have a strong hnRNPA1 site located within 10 nt (R _i  ≥ R _sequence ) [38] to ensure accurate proofreading of the exon [37].

Next, variants affecting the strength of SRFs were analyzed by a contextual exon definition analysis of ΔR _i,total . The context refers to the documented splicing activity of an SRF. For example, TIA1 has been shown to be an intronic enhancer of exon definition, so only intronic sites were considered. Similarly, hnRNPA1 proofreads the 3’ SS (acceptor) and inhibits exon recognition elsewhere [95]. Variants that lead to redundant SRFBS changes (i.e. one site is abolished and another proximate site [≤2 nt] of equivalent strength is activated) were assumed to have a neutral effect on splicing. If the strength of a site bound by PTB (polypyrimidine tract binding protein) was affected, its impact on binding by other factors was analyzed, as PTB impedes binding of other factors with overlapping recognition sites, but does not directly enhance or inhibit splicing itself [96].

To determine effects of variants on transcription factor (TF) binding, we first established which TFs bound to the sequenced regions of the gene promoters (and first exons) in this study by using ChIP-seq data from 125 cell types (Additional file 1: Methods) [97]. We identified 141 TFs with evidence for binding to the promoters of the genes we sequenced, including c-Myc, C/EBPβ, and Sp1, shown to transcriptionally regulate BRCA1, TP53, and ATM, respectively [98–100]. Furthermore, polymorphisms in TCF7L2, known to bind enhancer regions of a wide variety of genes in a tissue-specific manner [101], have been shown to increase risk of sporadic [102] and hereditary breast [103], as well as other types of cancer [104, 105].

IT-based models of the 141 TFs of interest were derived by entropy minimization of the DNase accessible ChIP-seq subsets [106]. Details are provided in Lu R, Mucaki E, and Rogan PK (BioRxiv; http://dx.doi.org/10.1101/042853). While some data sets would only yield noise or co-factor motifs (i.e. co-factors that bind via tethering, or histone modifying proteins [107]), techniques such as motif masking and increasing the number of Monte Carlo cycles yielded models for 83 TFs resembling each factor’s published motif. Additional file 3: Table S1 contains the final list of TFs and the models we built (described below) [108].

These TFBS models (N = 83) were used to scan all variants called in the promoter regions (10 kb upstream of transcriptional start site to the end of IVS1) of HBOC genes for changes in R _i . Binding site changes that weaken interactions with the corresponding TF (to R _i  ≤ R _sequence ) are likely to affect regulation of the adjacent target gene. Stringent criteria were used to prioritize the most likely variants and thus only changes to strong TFBSs (R _i,initial  ≥ R _sequence ), where reduction in strength was significant (ΔR _i  ≥ 4.0 bits), were considered. Alternatively, novel or strengthened TFBSs were also considered sources of dysregulated transcription. These sites were defined as having R _i,final  ≥ R _sequence and as being the strongest predicted site in the corresponding genomic interval (i.e. exceeding the R _i values of adjacent sites unaltered by the variant). Variants were not prioritized if the TF was known to a) enhance transcription and IT analysis predicted stronger binding, or b) repress transcription and IT analysis predicted weaker binding.

Two complementary strategies were used to assess the possible impact of variants within UTRs. First, SNPfold software was used to assess the effect of a variant on 2° structure of the UTR (Additional file 1: Methods) [20]. Variants flagged by SNPfold with the highest probability of altering stable 2° structures in mRNA (where p-value < 0.1) were prioritized. To evaluate these predictions, oligonucleotides containing complete wild-type and variant UTR sequences (Additional file 4: Table S2) were transcribed in vitro and followed by SHAPE analysis, a method that can confirm structural changes in mRNA [44].

Second, the effects of variants on the strength of RBBSs were predicted. Frequency-based, position weight matrices (PWMs) for 156 RNA-binding proteins (RBPs) were obtained from the RNA-Binding Protein DataBase (RBPDB) [109] and the Catalog of Inferred Sequence Binding Preferences of RNA binding proteins (CISBP-RNA) [110, 111]. These were used to compute information weight matrices (based on the method described by Schneider et al. 1984; N = 147) (see Additional file 1: Methods) [40]. All UTR variants were assessed using a modified version of the Shannon Pipeline [91] containing the RBPDB and CISBP-RNA models. Results were filtered to include a) variants with |ΔR _i | ≥ 4.0 bits, b) variants creating or strengthening sites (R _i,final  ≥ R _sequence and the R _i , _initial  < R _sequence ), and c) RBBSs not overlapping or occurring within 10 nt of a stronger, pre-existing site of another RBP.

Exonic protein-altering variant analysis

The predicted effects of all coding variants were assessed with SNPnexus [112–114], an annotation tool that can be applied to known and novel variants using up-to-date dbSNP and UCSC human genome annotations. Variants predicted to cause premature protein truncation were given higher priority than those resulting in missense (or synonymous) coding changes. Missense variants were first cross referenced with dbSNP 142 [115]. Population frequencies from the Exome Variant Server [116] and 1000Genomes [117] are also provided. The predicted effects on protein conservation and function of the remaining variants were evaluated by in silico tools: PolyPhen-2 [118], Mutation Assessor (release 2) [119, 120], and PROVEAN (v1.1.3) [121, 122]. Default settings were applied and in the case of PROVEAN, the “PROVEAN Human Genome Variants Tool” was used, which includes SIFT predictions as a part of its output. Variants predicted by all four programs to be benign were less likely to have a deleterious impact on protein activity; however this did not exclude them from mRNA splicing analysis (described above in IT-Based Variant Analysis). All rare and novel variants were cross-referenced with general mutation databases (ClinVar [123, 124], Human Gene Mutation Database [HGMD] [125, 126], Leiden Open Variant Database [LOVD] [127–134], Domain Mapping of Disease Mutations [DM²] [135], Expert Protein Analysis System [ExPASy] [136] and UniProt [137, 138]), and gene-specific databases (BRCA1/2: the Breast Cancer Information Core database [BIC] [139] and Evidence-based Network for the Interpretation of Germline Mutant Alleles [ENIGMA] [140]; TP53: International Agency for Research on Cancer [IARC] [141]), as well as published reports to prioritize them for further workup.

Variant classification

Flagged variants were prioritized if they were likely to encode a dysfunctional protein (indels, nonsense codon > 50 amino acids from the C-terminus, or abolition of a natural SS resulting in out-of-frame exon skipping) or if they exceeded established thresholds for fold changes in binding affinity based on IT (see Methods above). In several instances, our classification was superseded by previous functional or pedigree analyses (reported in published literature or databases) that categorized these variants as pathogenic or benign.

Positive control

Variant validation

Protein-truncating, prioritized splicing, and selected prioritized missense variants were verified by Sanger sequencing. Primers of PCR amplicons are indicated in Additional file 5: Table S3.

Deletion analysis

Capture, sequencing, and alignment

The average coverage of capture region per individual was 90.8x (range of 53.8 to 118.2x between 32 samples) with 98.8 % of the probe-covered nucleotides having ≥ 10 reads. Samples with fewer than 10 reads per nucleotide were re-sequenced and the results of both runs were combined. The combined coverage of these samples was, on average, 48.2x (±36.2).

The consistency of both library preparation and capture protocols was improved from initial runs, which significantly impacted sequence coverage (Additional file 1: Methods). Of the 102 patients tested, 14 had been previously Sanger sequenced for BRCA1 and BRCA2 exons. Confirmation of previously discovered SNVs served to assess the methodological improvements introduced during NGS and ultimately, to increase confidence in variant calling. Initially, only 15 of 49 SNVs in 3 samples were detected. The detection rate of SNVs was improved to 100 % as the protocol progressed. All known SNVs (N = 157) were called in subsequent sequencing runs where purification steps were replaced with solid phase reversible immobilization beads and where RNA bait was transcribed the same day as capture. To minimize false positive variant calls, sequence read data were aligned with CASAVA and CRAC, variants were called for each alignment with GATK, and discrepancies were then resolved by manual review.

GATK called 14,164 unique SNVs and 1147 indels. Only 3777 (15.3 %) SNVs were present in both CASAVA and CRAC-alignments for at least one patient, and even fewer indel calls were concordant between both methods (N = 110; 6.2 %). For all other SNVs and indels, CASAVA called 6871 and 1566, respectively, whereas CRAC called 13,958 and 110, respectively. Some variants were counted more than once if they were called by different alignment programs in two or more patients. Intronic and intergenic variants proximate to low complexity sequences tend to generate false positive variants due to ambiguous alignment, a well known technical issue in short read sequence analysis [146, 147], contributing to this discrepancy. For example, CRAC correctly called a 19 nt deletion of BRCA1 (rs80359876; also confirmed by Sanger sequencing) but CASAVA flagged the deleted segment as a series of false-positives (Additional file 6: Figure S1). For these reasons, all variants were manually reviewed.

IT-based variant identification and prioritization

Natural SS variants

The Shannon Pipeline reported 99 unique variants in natural donor or acceptor SSs. After technical and frequency filtering criteria were applied, 12 variants remained (Additional file 7: Table S4). IT analysis allowed for the prioritization of 3 variants, summarized in Table 2.

Patient ID	Gene	mRNA	rsID (dbSNP 142)	Information Change			Consequencef or Binding Factor Affected
				R i,initial	R i,final	ΔR i or R i e
			Allele Frequency (%)d	(bits)	(bits)	(bits)
Abolished Natural SS
7-4 F	ATM	c.3747-1G > Aa	Novel	11.0	0.1	−10.9	Exon skipping and use of alternative splice forms
4-1 F	ATM	c.6347 + 1G > Tb	Novel	10.4	−8.3	−18.6	Exon skipping
Leaky Natural SS
4-2B	CHEK2	c.320-5 T > Aa	rs121908700	6.8	4.1	−2.7	Leaky splicing with intron inclusion
			0.08
Activated Cryptic SS
7-3E	BRCA1	c.548-293G > A	rs117281398	−12.1	2.6	14.7	Cryptic site not expected to be used. Total information for natural exon is stronger than cryptic exon.
			0.74
7-4A	BRCA2	c.7618-269_7618-260del10	Novel	3.9	9.4	5.5	Cryptic site not expected to be used. Total information for natural exon is stronger than cryptic exon.
Pseudoexon formation due to activated acceptor SS
7-3 F	BRCA2	c.8332-805G > A	Novel	−9.3	5.4	5.6e	6065/211/592f
7-3D	CDH1	c.164-2023A > G	rs184740925	−6.6	4.3	6.5e	61,236/224/1798f
			0.3
5-3H	CDH1	c.2296-174 T > A	rs565488866	7.3	8.5	5.0e	1175/50/124f
			0.02
Pseudoexon formation due to activated donor SS
3-6A	BRCA1	c.212 + 253G > A	rs189352191	4.1	6.7	5.2e	186/63/1250f
			0.08
5-2G	BRCA2	c.7007 + 2691G > A	rs367890577	4.7	7.2	7.7e	2589/103/5272f
			0.02
Affected TFBSs
7-4B	BRCA1	c.-8895G > A	Novel	10.9	−0.2	−11.1	GATA-3 (GATA3)
5-3E	CDH1	c.-54G > C	rs5030874	1.7	12.0	10.4	E2F-4 (E2F4)
7-4E			0.16
5-2B	PALB2	c.-291C > G	rs552824227	12.1	−1.3	−13.4	GABPα (GABPA)
			0.1
7-2 F	TP53	c.-28-3132 T > C	rs17882863	−6.3	10.9	17.2	RUNX3 (RUNX3)
			0.3
4-1A	TP53	c.-28-1102 T > C	rs113451673	5.1	12.3	7.2	E2F-4 (E2F4)
			0.4	8.0	12.9	4.8	Sp1 (SP1)
Affected RBBSs
7-4G	ATM	c.-244 T > A	rs539948218	9.8	−19.9	−29.7	RBFOX
		c.-744 T > A	0.04
		c.-1929 T > A
		c.-3515 T > A
5-3C	CDH1	c.*424 T > A	Novel	−20.3	9.6	29.9	SF3B4
				8.2	1.8	−6.4	CELF4
7-2E	CHEK2	c.-588G > A	rs141568342	10.9	3.7	−7.2	BX511012.1
4-3C.5-4G	CHEK2	c.-345C > Tc	rs137853007	3.3	11.4	8.2	SF3B4
3-1A	TP53	c.-107 T > C	rs113530090	10.5	4.5	−6.0	ELAVL1
4-1H		c.-188 T > C	0.72
4-2H	TP53	c.*1175A > C	rs78378222	10.7	4.1	−6.6	KHDRBS1
7-2 F		c.*1376A > C	0.26
		c.*1464A > C

^aConfirmed by Sanger sequencing

^bAmbiguous Sanger sequencing results

^cPrioritized under missense change and was therefore verified with Sanger sequencing. Variant was confirmed

^dIf available

^e R _i of site of opposite polarity in the pseudoexon

^fConsequences for pseudoexon formation describe how the intron is divided: “new intron A length/pseudoexon length/new exon B length

None of the variants have been previously reported by other groups with the exception of CHEK2 c.320-5T>A [148]

First, the novel ATM variant c.3747-1G > A (chr11:108,154,953G > A; sample number 7-4 F) abolishes the natural acceptor of exon 26 (11.0 to 0.1 bits). ASSEDA reports the presence of a 5.3 bit cryptic acceptor site 13 nt downstream of the natural site, but the effect of the variant on a pre-existing cryptic site is negligible (~0.1 bits). The cryptic exon would lead to exon deletion and frameshift (Fig. 3a). ASSEDA also predicts skipping of the 246 nt exon, as the R _i,final of the natural acceptor is now below R _i,minimum (1.6 bits), altering the reading frame. Second, the novel ATM c.6347 + 1G > T (chr11:108188249G > T; 4-1 F) abolishes the 10.4 bit natural donor site of exon 44 (ΔR _i  = -18.6 bits), and is predicted to cause exon skipping. Finally, the previously reported CHEK2 variant, c.320-5A > T (chr22:29,121,360 T > A; rs121908700; 4-2B) [148] weakens the natural acceptor of exon 3 (6.8 to 4.1 bits), and may activate a cryptic acceptor (7.4 bits) 92 nt upstream of the natural acceptor site which would shift the reading frame (Fig. 4). A constitutive, frameshifted alternative isoform of CHEK2 lacking exons 3 and 4 has been reported, but skipping of exon 3 alone is not normally observed.

Variants either strengthening (N = 4) or slightly weakening (ΔR _i  < 1.0 bits; N = 4) a natural site were not prioritized. In addition, we rejected the ATM variant (c.1066-6 T > G; chr11:108,119,654 T > G; 4-1E and 7-2B), which slightly weakens the natural acceptor of exon 9 (11.0 to 8.1 bits). Although other studies have shown leaky expression as a result of this variant [149], a more recent meta-analysis concluded that this variant is not associated with increased breast cancer risk [150].

UTR structure and protein binding

There were 364 unique UTR variants found by sequencing. These variants were evaluated for their effects on mRNA 2° structure (including that of splice forms with alternate UTRs in the cases of BRCA1 and TP53) through SNPfold, resulting in 5 flagged variants (Table 3), all of which have been previously reported.

Classa	Patient ID	Gene	mRNA	UTR position	rsID (dbSNP 142)	Ranke	p-value
					Allele Frequency (%)d
F	In 26 patients	BRCA2 b	c.-52A > G	5’ UTR	rs206118	2/900	0.002
					14.86
F	In 40 patients	BRCA2 b	c.*532A > G	3’ UTR	rs11571836	239/2700	0.089
					19.75
P	7-4C	CDH1 c	c.-71C > G	5’ UTR	rs34033771	69/600	0.115
					0.56
F	4-2E	TP53 b	c.*485G > A	3’ UTR	rs4968187	169/4500	0.038
	5-4A
					5.11
F	2-1A, 7-1B, 5-2A.7-1D, 7-2B, 7-2F	TP53 b	c.*826G > A	3’ UTR	rs17884306	371/4500	0.082
	7-4C
					5.71

^aF:Flagged; P:Prioritized

^bLong Range UTR SNPfold Analysis

^cLocal Range SNPfold Analysis

^dIf available

^eRank of the SNP, in terms of how much it changes the mRNA structure compared to all other possible mutations

Analysis of three variants using mFOLD [83] revealed likely changes to the UTR structure (Fig. 5). Two variants with possible 2° structure effects were common (BRCA2 c.-52A > G [N = 26 samples] and c.*532A > G [N = 40]) and not prioritized. The 5’ UTR CDH1 variant c.-71C > G (chr16:68771248C > G; rs34033771; 7-4C) disrupts a double-stranded hairpin region to create a larger loop structure, thus increasing binding accessibility (Fig. 5a and b). Analysis using RBPDB and CISBP-RNA-derived IT models suggests this variant affects binding by NCL (Nucleolin, a transcription coactivator) by decreasing binding affinity 14-fold (R _i,initial  = 6.6 bits, ΔR _i  = -3.8 bits) (Additional file 10: Table S7). This RBP has been shown to bind to the 5’ and 3’ UTR of p53 mRNA and plays a role in repressing its translation [151].

In addition, the TP53 variant c.*485G > A (NM_000546.5: chr17:7572442C > T; rs4968187) is found at the 3’ UTR and was identified in two patients (4-2E and 5-4A). In silico mRNA folding analysis demonstrated this variant disrupts a G/C bond of a loop in the highest ranked potential mRNA structure (Fig. 5c and d). Also, SHAPE analysis showed a difference in 2° structure between the wild-type and mutant (data not shown). IT analysis with RBBS models indicated that this variant significantly increases the binding affinity of SF3B4 by > 48-fold (R _i,final  = 11.0 bits, ΔR _i  = 5.6 bits) (Additional file 10: Table S7). This RBP is one of four subunits comprising the splice factor 3B, which binds upstream of the branch-point sequence in pre-mRNA [152].

The third flagged variant also occurs in the 3’ UTR of TP53 (c.*826G > A; chr17:7572,101C > T; rs17884306), and was identified in 6 patients (2-1A, 7-1B, 5-2A.7-1D, 7-2B, 7-2 F, and 7-4C). It disrupts a potential loop structure, stabilizing a double-stranded hairpin, and possibly making it less accessible (Fig. 5e and f). Analysis using RBPDB-derived models suggests this variant could affect the binding of both RBFOX2 and SF3B4 (Additional file 10: Table S7). A binding site for RBFOX2, which acts as a promoter of alternative splicing by favoring the inclusion of alternative exons [153], is created (R _i,final  = 9.8 bits; ΔR _i  = -6.5 bits). This variant is also expected to simultaneously abolish a SF3B4 binding site (R _i,final  = -20.3 bits; ΔR _i  = -29.9 bits).

RBPDB- and CISBP-RNA-derived information model analysis of all UTR variants resulted in the prioritization of 1 novel, and 5 previously-reported variants (Table 2). No patient within the cohort exhibited more than one prioritized RBBS variant.

To evaluate the background rate of prioritizing variants flagged by this method, all 5’ and 3’ UTR SNVs in dbSNP144 for the 7 genes sequenced (excluding those already flagged in Table 3) were evaluated by SNPfold and our RBP information models. Of 1207 SNVs, only 10 were prioritized with both methods, which results in a background rate of 0.83 %.

Exonic variants altering protein sequence

Exonic variants called by GATK (N = 245) included insertions, deletions, nonsense, missense, and synonymous changes.

Variant classification

Initially, 15,311 unique variants were identified by complete gene sequencing of 7 HBOC genes. Of these, 132 were flagged after filtering, and further reduced by IT-based variant analysis and consultation of the published literature to 87 prioritized variants. Figure 6 illustrates the decrease in the number of unique variants per patient at each step of our identification and prioritization process. The distribution of prioritized variants by gene is 34 in ATM, 13 in BRCA1, 11 in BRCA2, 8 in CDH1, 6 in CHEK2, 10 in PALB2, and 5 in TP53 (Additional file 13: Table S10), which are categorized by type in Table 5.

	Indel	Nonsense	Missense	Natural Splicing	Cryptic Splicing	Pseudoexon	SR Factor	TF	UTR Structure	UTR Binding	Total
ATM	0	0	14	2	0	0	18	0	0	1	34 a
BRCA1	2	0	2	0	0	1	7	1	0	0	13
BRCA2	0	2	3	0	0	2	4	0	0	0	11
CDH1	0	0	2	0	0	2	1	1	1	1	8
CHEK2	0	0	2	1	0	0	3	0	0	2	6 a
PALB2	1	2	3	0	0	0	3	1	0	0	10
TP53	0	0	1	0	0	0	0	2	0	2	5
Total	3	4	27	3	0	5	36	5	1	6

Three variants were prioritized under multiple categories: ATM chr11:108121730A > G (missense and SRFBS), CHEK2 chr22:29121242G > A (missense, UTR binding), and CHEK2 chr22:29130520C > T (missense, UTR binding)

^a Counts represent the number of unique variants identified (i.e. a variant is not counted twice if it appeared in multiple individuals)

Three prioritized variants have multiple predicted roles: ATM c.1538A > G in missense and SRFBS, CHEK2 c.190G > A in missense and UTR binding, and CHEK2 c.433C > T in missense and UTR binding. Of the 102 patients that were sequenced, 72 (70.6 %) exhibited at least one prioritized variant, and some patients harbored more than one prioritized variant (N = 33; 32 %). Additional file 14: Table S11 presents a summary of all flagged and prioritized variants for patients with at least one prioritized variant.

Prioritization of potential deletions

Using BreakDancer, none of the individuals analyzed exhibited large rearrangements that met the level of stringency required, but a small intragenic rearrangement in BRCA1 was identified and confirmed by Sanger sequencing. Attempts to detect deletions with BreakDancer only flagged single, non-contiguous paired-end reads, rather than a series of reads clustered within the same region within the same individual, which would be necessary to indicate the presence of a true deletion or structural rearrangement.

After prioritizing individuals for potential hemizygosity in the sequenced regions, potential deletions were detected in BRCA2 and CDH1. Patient UWO5-4D exhibited a non-polymorphic 32.1 kb interval in BRCA2, spanning introns 1 to 13, that was absent from all of the other individuals (chr13:32890227-32922331). Haploview (hapmap.org) showed very low levels of LD in this region. The potential deletion may extend further downstream, however the presence of a haploblock, covering the entire sequenced interval beyond exon 11, with significant LD precludes delineation of the telomeric breakpoint. We also flagged a non-polymorphic 2.6 kb interval near the 3’ end of CDH1 in 6 individuals (UWO3-5, UWO4-2C, UWO4-4E, UWO4-4 F, UWO4-2G, UWO5-2H). This is a low LD region spanning chr16:68861286-68863887 that includes exons 14 and 15, and is polymorphic in all of the other individuals sequenced. CDH1 mutations are characteristically present in families with predisposition to gastric cancer, however breast cancer frequently co-occurs [69]. A study of CDH1 deletions in inherited gastric cancer identified two families with deletions that overlap the intervals prioritized in the present study [162].

Comparison to combined annotation dependent depletion

The analysis and prioritization of non-coding variants can also be accomplished using Combined Annotation Dependent Depletion (CADD; [163]), which uses known and simulated variants to compute a C-score, an ad hoc measure of how deleterious is likely to be. The suggested C-score cutoff is between 10 and 20, though it is stated that any selected cutoff value would be arbitrary (http://cadd.gs.washington.edu/info). This contrasts with information-based methods, which are based on thermodynamically-defined thresholds. To directly compare methods, CADD scores were obtained for all prioritized or flagged SNVs. Half of prioritized variants met this cutoff (C > 10), while only 28.6 % of flagged variants did the same. All prioritized nonsense variants (4/4) and 26/27 missense variants had strong C-scores. Prioritized non-coding variant categories that correlated well with CADD include natural splicing variants (4/4), UTR structure variants (1/1), and RBPs (4/6). Weakly correlated variants included those affecting SRFBPs (5/36), TFBS (2/5), and pseudoexon activating variants (0/5). Missense mutations comprised 75 % of the flagged variants with C > 10. The aforementioned flagged splicing variant ATM c.1066-6 T > G also exceeded the threshold C value (C = 11.9). Meanwhile, the flagged TP53 variant, shown by SHAPE analysis to alter UTR structure, did not (C = 5.3). Despite consistency between some variant categories, the underlying assumptions of each approach probably explain why these results differ for non-coding variants. The limited numbers validated, deleterious non-coding variants also contributes to the accuracy of these predictions [163].

Variant verification

We verified prioritized protein-truncating (N = 7) and splicing (N = 4) variants by Sanger sequencing (Tables 2 and 4, respectively). In addition, two missense variants (BRCA2 c.7958 T > C and CHEK2 c.433C > T) were re-sequenced, since they are indicated as likely pathogenic/pathogenic in ClinVar (Additional file 12: Table S9). All protein-truncating variants were confirmed, with one exception (BRCA2 c.7558C > T, no evidence for the variant was present for either strand). Two of the mRNA splicing mutations were confirmed on both strands, while the other two were confirmed on a single strand (ATM c.6347 + 1G > T and ATM c.1066-6 T > G). Both documented pathogenic missense variants were also confirmed.

Non-coding variants

Although coding variants are typically the sole focus of a molecular diagnostic laboratory (with the exception of the canonical dinucleotide positions within SS), non-coding mutations have long been known to be disease causing [19, 36, 175–183]. In this study, variant density in non-coding regions significantly exceeded exonic variants by > 60-fold, which, in absolute terms, constituted 1.6 % of the 15,311 variants. This is comparable to whole genome sequencing studies, which typically result in 3-4 million variants per individual, with < 2 % occurring in protein coding regions [184]. IT analysis prioritized 3 natural SS, 36 SRFBS, 5 TFBS, and 6 RBBS variants and 5 predicted to create pseudoexons. Two SS variants in ATM (c.3747-1G > A and c.6347 + 1G > T) were predicted to completely abolish the natural site and cause exon skipping. A CHEK2 variant (c.320-5A > T) was predicted to result in leaky splicing.

The IT-based framework evaluates all variants on a common scale, based on bit values, the universal unit that predicts changes in binding affinity [185]. A variant can alter the strength of one or a “set” of binding sites; the magnitude and direction of these changes is used to rank their significance. The models used to derive information weight matrices take into account the frequency of all observed bases at a given position of a binding motif, making them more accurate than consensus sequence and conservation-based approaches [36].

IT has been widely used to analyze natural and cryptic SSs [36], but its use in SRFBS analysis was only introduced recently [38]. For this reason, we assigned conservative, minimum thresholds for reporting information changes. Although there are examples of disease-causing variants resulting in small changes in R _i [174, 186–192], the majority of deleterious splicing mutations that have been verified functionally, produce large information changes. Among 698 experimentally verified deleterious variants in 117 studies, only 1.96 % resulted in < 1.0 bit change [36]. For SRFBS variants, the absolute information changes for deleterious variants ranged from 0.2 to 17.1 bits (mean 4.7 ± 3.8). This first application of IT in TFBS and RBBS analysis, however, lacks a large reference set of validated mutations for the distribution of information changes associated with deleterious variants. The release of new ChIP-seq datasets will enable IT models to be derived for TFs currently unmodeled and will improve existing models [193].

Pseudoexon activation results in disease-causing mutations [194], however such consequences are not customarily screened for in mRNA splicing analysis. IT analysis was used to detect variants that predict pseudoexon formation and 5 variants were prioritized. Previously, we have predicted experimentally proven pseudoexons with IT (Ref 42: Table 2, No #2; and Ref 195: Table 2, No #7) [42, 195]. Although it was not possible to confirm prioritized variants in the current study predicted to activate pseudoexons because of their low allele frequencies, common intronic variants that were predicted to form pseudoexons were analyzed. We then searched for evidence of pseudoexon activation in mapped human EST and mRNA tracks [196] and RNA-seq data of breast normal and tumour tissue from the Cancer Genome Atlas project [15]. One of these variants (rs6005843) appeared to splice the human EST HY160109 [197] at the predicted cryptic SS and is expressed within the pseudoexon boundaries.

Variants that were common within our population sample (i.e. occurring in > 5 individuals) and/or common in the general population (>1.0 % allele frequency) reduced the list of flagged variants substantially. This is now a commonly accepted approach for reducing candidate disease variants [166], based on the principle that the disease-causing variants occur at lower population frequencies. Variants occurring in > 5 patients all either had allele frequencies above 1.0 % or, as shown previously, resulted in very small ΔR _i values [198].

The genomic context of sequence changes can influence the interpretation of a particular variant [36]. For example, variants causing significant information changes may be interpreted as inconsequential if they are functionally redundant or enhancing existing binding site function (see IT-Based Variant Analysis for details). Our understanding of the roles and context of these cognate protein factors is incomplete, which affects confidence in interpretation of variants that alter binding. Also, certain factors with important roles in the regulation of these genes, but that do not bind DNA directly or in a sequence-specific manner (eg. CtBP2 [199]), could not be included. Therefore, some variants may have been incorrectly excluded.

Prioritization of potential deletions

Although individuals can be prioritized based on potential hemizygosity, this does not definitively identify deletions. Nevertheless, it should be possible to prioritize those individuals worthy of further detailed diagnostic workup. It has not escaped our attention that the weighted probabilities obtained from this analysis could be represented and formalized using the same units of Shannon information (in bits) as the other sequence changes we have described, analogous to single or multinucleotide gene variants predicted to affect nucleic acid binding sites. Full development and validation of this method is in progress.

Coding sequence changes

We also identified 4 nonsense and 3 indels in this cohort. In one individual, a 19 nt BRCA1 deletion in exon 15 causes a frameshift leading to a stop codon within 14 codons downstream. This variant, rs80359876, is considered clinically relevant. Interestingly, this deletion overlaps two other published deletions in this exon (rs397509209 and rs80359884). This raises the question as to whether this region of the BRCA1 gene is a hotspot for replication errors. DNA folding analysis indicates a possible 15 nt long stem-loop spanning this interval as the most stable predicted structure (data not shown). This 15 nt structure occurs entirely within the rs80359876 and rs397509209 deletions and partially overlaps rs80359884 (13 of 15 nt of the stem loop). It is plausible that the 2° structure of this sequence predisposes to a replication error that leads to the observed deletion.

Missense coding variants were also assessed using multiple in silico tools and evaluated based on allele frequency, literature references, and gene-specific databases. Of the 27 prioritized missense variants, the previously reported CHEK2 variant c.433G > A (chr22:29121242G > A; rs137853007) stood out, as it was identified in one patient (4-3C.5-4G) and is predicted by all 4 in silico tools to have a damaging effect on protein function. Accordingly, Wu et al. (2001) demonstrated reduced in vitro kinase activity and phosphorylation by ATM kinase compared to the wild-type CHEK2 protein [200], presumably due to the variant’s occurrence within the forkhead homology-associated domain, involved in protein-phosphoprotein interactions [201]. Implicated in Li-Fraumeni syndrome, known to increase the risk of developing several types of cancer including breast [202, 203], the CHEK2: c.433G > A variant is expected to result in a misfolded protein that would be targeted for degradation via the ubiquitin-proteosome pathway [204]. Another important missense variant is c.7958 T > C (chr13:32,936,812 T > C; rs80359022; 4-4C) in exon 17 of BRCA2. Although classified as being of unknown clinical importance in both BIC and ClinVar, it has been classified as pathogenic based on posterior probability calculations [205].

It is unlikely that all prioritized variants are pathogenic in patients carrying more than one prioritized variant. Nevertheless, a polygenic model for breast cancer susceptibility, whereby multiple moderate and low-risk alleles contribute to increased risk of HBOC may also account for multiple prioritized variants [206, 207]. There was a significant fraction of patients (29.4 %) in whom no variants were prioritized. This could be due to a) the inability of the analysis to predict a variant affecting the binding sites analyzed, b) a pathogenic variant affecting a function that was not analyzed or in a gene that was not sequenced, c) a large rearrangement/deletion where both breakpoints occur beyond the captured genomic intervals (which is unlikely, as this would have been observed as an extended non-polymorphic sequence), or d) the significant family history was not due to heritable, but instead to shared environmental influences.

BRCA coding variants were found in individuals who were previously screened for lesions in these genes, suggesting this NGS protocol is a more sensitive approach for detecting coding changes. However, previous testing of a number of these patients had been predominantly based on PTT and MLPA, which have lower sensitivity for detecting mutations than sequence analysis. Nevertheless, we identified 2 BRCA1 and 2 BRCA2 variants predicted to encode prematurely truncated proteins. Fewer non-coding BRCA variants were prioritized (15.7 %) than expected by linkage analysis [49], however this presumes at least 4 affected breast cancer diagnoses per pedigree, and, in the present study, the number of affected individuals per family was not known.

Prioritization of a variant does not equate with pathogenicity. Some prioritized variants may not increase risk, but may simply modify a primary unrecognized pathogenic mutation. A patient with a known BRCA1 nonsense variant, used as a positive control, was also found to possess an additional prioritized variant in BRCA2 (missense variant chr13:32911710A > G), which was flagged by PROVEAN and SIFT as damaging, as well as flagged for changing an SRFBS for abolishing a PTB site (while simultaneously abolished an exonic hnRNPA1 site). This variant has been identified in cases of early onset prostate cancer and is considered a VUS in ClinVar [143]. Similarly, variants prioritized in multiple patients may act as risk modifiers rather than pathogenic mutations. A larger cohort of patients with known pathogenic mutations would be necessary to calculate a background/basal rate of falsely flagged variants.

Other groups have attempted to develop comprehensive approaches for variant analysis, analogous to the one proposed here [208–210]. While most employ high-throughput sequencing and classify variants, either the sequences analyzed or the types of variants assessed tend to be limited. In particular, non-coding sequences have not been sequenced or studied to the same extent, and none of these analytical approaches have adopted a common framework for mutation analysis.

Our published oligonucleotide design method [77] produced an average sequence coverage of 98.8 %. The capture reagent did not overlap conserved highly repetitive regions, but included divergent repetitive sequences. Nevertheless, neighboring probes generated reads with partial overlap of repetitive intervals. As previously reported [147], we noted that false positive variant calls within intronic and intergenic regions were the most common consequence of dephasing in low complexity, pyrimidine-enriched intervals. This was not alleviated by processing data with software programs based on different alignment or calling algorithms. Manual review of all intronic or intergenic variants became imperative. As these sequences can still affect functional binding elements detectable by IT analysis (i.e. 3’ SSs and SRFBSs), it may prove essential to adopt or develop alignment software that explicitly and correctly identifies variants in these regions [147]. Most variants were confirmed with Sanger sequencing (10/13), and those that could not be confirmed are not necessarily false positives. A recent study demonstrated that NGS can identify variants that Sanger sequencing cannot, and reproducing sequencing results by NGS may be worthwhile before eliminating such variants [211].

Availability of supporting data

Variants will be deposited with the ENIGMA Consortium (www.enigmaconsortium.org), which is a designated organization for curation of HBOC mutations and which is charged with protection of genetic privacy of participants. Additionally, all likely pathogenic variants were submitted to ClinVar (submission ID: SUB1332053) while other novel variants were submitted to dbSNP (ss1966658584-1966658622).