This study was approved and carried out in accordance with the recommendations of the Western Institutional Review Board (WIRB). All subjects gave written informed consent in accordance with the Declaration of Helsinki. Blood from normal healthy controls or patients participating in the Ambry Genetics Family Studies program was drawn in PAXgene Blood RNA Tubes and stored according to the manufacturer's recommendations (PreAnalytiX, Hombrechtikon, Switzerland). RNA was extracted using the PAXgene Blood RNA Kit according to the recommended protocol (PreAnalytiX). Informed consent was obtained from all participants. Breast RNA was purchased from Amsbio (Lake Forest, CA, USA) and BioChain (Newark, CA, USA). RNA quality was determined using the RIN number calculated by the TapeStation 2200 with RNA ScreenTape or High Sensitivity RNA ScreenTape (Agilent, Santa Clara, CA, USA).

Lymphoblastoid cell lines (LCLs) were obtained from the Kathleen Cuningham Consortium for Research into Familial Breast Cancer (kConFab, Melbourne, Australia) from 4 carriers of BRCA1 or BRCA2 variants and 2 controls. Genotypes were verified by Sanger sequencing. LCLs were maintained according to the recommendations of kConFab. Inhibition of nonsense-mediated decay (NMD) was performed using puromycin (300 μg/ml) or cycloheximide (100 μg/ml) for 4 h, as previously described (8, 14).

cDNA was generated using the SuperScript IV First-Strand Synthesis System (Thermo Fisher Scientific, Chino, CA, USA). PCR was performed using either Platinum SuperFi PCR Master Mix (Thermo Fisher Scientific) or HotStarTaq Master Mix (Qiagen, Valencia, CA, USA) as previously described (8).

PCR products were analyzed using digital electrophoresis with D1000 ScreenTape and Reagents on the TapeStation 2200 (Agilent). Capillary electrophoresis (CE) was performed on an ABI 3730xl using MapMarker1000 as a standard (BioVentures, Murfreesboro, TN, USA). Primers were tagged at the 5′ end with FAM or HEX for detection by CE. CE analysis was performed with GeneMapper software (Thermo Fisher Scientific).

PCR products were cloned into pGEM-T Easy and transformed into bacteria according to the manufacturer's recommended protocol (Promega, Fitchburg, WI, USA). Individual white colonies were picked, amplified by rolling-circle replication, and Sanger sequenced by Genewiz (La Jolla, CA, USA).

For CloneSeq, cDNA, PCR, and cloning were performed as described above. All colonies on a plate were scraped and suspended in PBS. Plasmids were extracted with the GeneJET Plasmid Miniprep kit (Thermo Fisher Scientific). CloneSeq libraries were constructed according to the protocol outlined by KAPA Biosystems (Wilmington, MA, USA) using the Hyper Prep kit. Briefly, DNA was sheared to an average size of 250–350 bp using sonication (Covaris, Woburn, MA, USA). DNA fragment ends were repaired and phosphorylated. An “A” base was added to the 3′ end of the blunted fragments, followed by ligation of single-indexed adapters via T-A mediated ligation. The size and concentration of the DNA library were determined using the TapeStation 2200. Massively Parallel Sequencing (MPS) was performed on an Illumina MiSeq, which generated 2 × 250 paired-end reads. Sequencing reads were aligned to the hg19 reference genome and analyzed using Ambry's Bioinformatic Pipeline (see below).

For whole transcriptome RNA-Seq, globin mRNA and ribosomal RNA were depleted using the Globin-Zero Gold rRNA removal kit (Illumina, San Diego, CA, USA). After depletion, RNA was fragmented and single-indexed cDNA libraries were generated using an RNA Hyper Prep kit (KAPA Biosystems). Quality control was performed using the TapeStation 2200. Libraries were checked for average fragment size, concentration, and the presence of spurious peaks such as adapter dimers. Concentration was confirmed using a Qubit Fluorometer (ThermoFisher). Libraries were sequenced to a depth of 1 × 10⁸ paired end reads (2 × 150 bp) per sample on the Illumina NextSeq platform. Sequencing reads were aligned to the hg19 reference genome and analyzed using Ambry's Bioinformatic Pipeline for alternative splicing events and differentially expressed genes.

Paired-end RNA-seq reads (2 × 250 bp) and Sanger sequencing reads (~1,100 bp) were first aligned to the hg19 human reference genome. For Sanger reads, GMAP aligner (version 2016-04-04) was used with default parameters to perform single transcript alignment (STA) of very long reads. For CloneSeq reads, STAR aligner v2.5.2a was used with default parameters except the “outSAMtype” parameter was set to “BAM SortedByCoordinate.” The mapped reads were then analyzed by our customized Ambry Bioinformatics Pipeline (ABP) software to detect splicing events such as exon skipping, alternative 5′ donor site, alternative 3′ acceptor site, and intron retention (15). These events are detected by the pipeline, based on the alignments against the reference genome (Figure 1B): (1) exon skipping, if there is no reads align to one exon or several consecutive exons; (2) partial exon skipping, if there is no read alignment in one end of an exon; (3) partial intron inclusion, if there is alignment in one end of an intron; (4) intron retention, if there is alignment in a whole intron; (5) cryptic exon, if there is alignment in the middle of an intron and no alignment in the rest of the intron. Schematic representations of these splicing events are illustrated in Figure 1B. To quantify splicing events, we calculated the percentage of alternative splicing event against a given transcript/isoform: percent of alternative splicing event = (number of reads supporting alternative splicing event)/(number of all reads in the region covering alternative splicing event). To filter out noise caused by sequencing and alignment errors, or due to the expression of ultra-rare isoforms, the splicing events with “number of reads supporting alternative splicing event” <20, or “number of all reads in the region covering splicing event” <50, or “percent of splicing event” <2.5% were filtered out. HGVS nomenclature values were approximate for intron retention and alternative splicing site events due to differences in alignments based on NGS reads.

The variant BRCA1 c.5467+5G>T, which impairs the native donor splice site of BRCA1 exon 23, has been described to result in skipping of exon 23 (Δ23) (8). This variant is currently classified as VUS by ENIGMA (class 3). LCLs were obtained from carriers of the variant and 2 controls, and reverse-transcriptase PCR (RT-PCR) was performed using the conditions recommended by Whiley et al. including the use of the same primers, reverse transcriptase, and NMD inhibitors puromycin (puro) and cycloheximide (CHX) (8). Skipping of exon 23 (Δ23) was clearly detected by digital electrophoresis of the RT-PCR products in the BRCA1 c.5467+5G>T carrier's LCL, but not in two control LCLs treated with NMD inhibitors (Figure 2A). RT-PCR products were cloned and CloneSeq was performed on these samples for sequence characterization and quantification. The sequence and absolute number of reads observed in the carrier and control cell lines are shown using Sashimi plots (Figure 2B). Sashimi plots provide a quantitative visualization of aligned MPS reads that enables quantitative comparison of exon usage across samples (16). A total of 4,018 reads supporting exon 23 skipping (r.5407_5467del61) were detected in the BRCA1 c.5467+5G>T LCL, whereas none was detected in the control LCL (Figure 2B). Abnormal transcripts levels were then measured as a “percent spliced in index” (PSI) (Figure 2C). PSI demonstrates the ratio between reads including or excluding exons, indicating how efficiently sequences of interest are spliced into transcripts (17). This analysis indicated that ~25% of transcripts expressed by the BRCA1 c.5467+5G>T LCL contains skipping of exon 23, whereas skipping of exon 23 was not detected in negative control LCLs (Figure 2C).

To validate the CloneSeq results we performed, in the same set of samples, the mRNA splicing assays recommended by the members of the ENIGMA consortium (8), including capillary electrophoresis (Figure 2D), and Sanger sequencing of subcloned transcripts (Figure 2E). Capillary electrophoresis clearly detected Δ23 in the BRCA1 c.5467+5G>T LCL, in addition to the full-length WT transcript (Figure 2D). Sanger sequencing also detected Δ23 exclusively in the BRCA1 c.5467+5G>T carrier's LCL (Figure 2E).

The variant BRCA1 c.135-1G>T, which impairs the native acceptor site of BRCA1 exon 5, has been associated with multiple splicing isoforms (8, 18), including an abundant transcript with skipping of exon 5 (Δ5). This variant is currently classified as pathogenic by ENIGMA (class 5). RT-PCR for the BRCA1 c.135-1G>T carrier's LCL and control LCLs was performed following ENIGMA recommendations and analyzed by digital electrophoresis, which detected a band consistent in size with Δ5 (Figure 3A). RT-PCR products were cloned and CloneSeq performed, and the sequence and absolute number of reads observed in the carrier and control cell lines are shown using Sashimi plots (Figure 3B). A total of 10,902 reads supporting exon 5 skipping (r.135_212del78) were detected in the BRCA1 c.135-1G>T LCL, whereas only 584 reads were detected in the control LCL (Figure 3B). Quantification of splicing events indicated that the BRCA1 c.135-1G>T LCL expresses ~50% of transcripts with skipping of exon 5, whereas LCL negative controls have negligible levels of Δ5 (Figure 3C). Individual colonies were selected for transcript confirmation by Sanger sequencing, and the Δ5 transcript was the most abundant abnormal transcripts in the BRCA1 c.135-1G>T LCL (Figure 3D).

BRCA2 c.8632+1G>A, located at intron 20 of BRCA2, was shown to result in skipping of exon 20 (Δ20) (8), and is currently classified by ENIGMA as a pathogenic alteration (class 5). RT-PCR of the BRCA2 c.8632+1G>A carrier's LCL and control LCLs was performed following ENIGMA recommendations. Digital electrophoresis analysis of the RT-PCR products detected several alternative transcripts, in addition to the full-length mRNA (Figure 4A). RT-PCR products were cloned and CloneSeq was performed. The sequence and absolute number of reads observed are shown in Sashimi plots for the carrier and control cell lines demonstrating a total of 3,002 reads supporting exon 20 skipping (r.8488_8632del145) in the BRCA2 c.8632+1G>A LCL, whereas only 25 reads supporting Δ20 were detected in the control LCL (Figure 4B). Quantification of splicing events indicated that ~20% of the transcripts expressed by the BRCA2 c.8632+1G>A LCL have skipping of exon 20, whereas LCL negative controls have negligible levels of Δ20 (Figure 4C). Several alternative spliced transcripts, detected both in the BRCA2 c.8632+1G>A LCL and controls, were also identified (Figure 4C). Individual colonies were selected for transcript confirmation by Sanger sequencing, which identified the most abundant abnormal transcript Δ20 in the BRCA2 c.8632+1G>A LCL, in addition to confirming other minor alternative isoforms detected by CloneSeq (Figure 4D).

BRCA2 c.9501+3A>T is located in the native donor site of intron 25. This variant was reported to result in low levels of skipping of exon 25 (Δ25), and it is currently classified as benign (class 1) by ENIGMA (8). RT-PCR was performed on the carrier's LCL and on the control cells. Digital electrophoresis identified a minor band consistent with the size of Δ25 in the BRCA2 c.9501+3A>T LCL that was not detected in negative controls (Figure 5A). RT-PCR products were cloned, and individual colonies were selected for Sanger sequencing, which detected Δ25 (r.9257_9501del245) in the BRCA2 c.9501+3A>T LCL, in addition to low levels of other alternatively spliced transcripts in controls (Figure 5B). To more accurately quantify the relative abundance of Δ25, we performed CloneSeq. The assay detected a total of 2,883 reads supporting exon 25 skipping in the BRCA2 c.9501+3A>T LCL, whereas no reads supporting skipping of this exon were detected in the control LCL (Figure 5C). The reads supporting Δ25 were ~10% of the total splicing events detected in the BRCA2 c.9501+3A>T LCL (Figure 5D).

Subsequently, we compared CloneSeq LCL results with an analysis of RNA isolated from the blood cells of carriers of the BRCA2 c.9501+3A>T alteration. RT-PCR performed on RNA from individuals that are heterozygous for BRCA2 c.9501+3A>T (proband and mother) and negative controls (father, LCL-, normal breast RNA, and normal blood controls) detected Δ25 only in the positive samples (Figure 5E). CloneSeq was performed to quantify Δ25 in these samples, which detected 9,279 and 4,066 splicing events supporting Δ25 in the proband and mother's samples respectively, while none were detected in the negative controls (Figure 5F). Interestingly, the percentage of Δ25 transcripts varies among different carriers. For the proband, Δ25 represents ~20% of splicing events, while it accounted for ~10% in the mother's sample (Figure 5G), ~5% in the LCL+ BRCA2 c.9501+3A>T without NMD inhibition (Figure 5G), and ~10% of splicing events when the LCL+ is treated with inhibitors (Figure 5D).

We next analyzed a novel uncharacterized VUS, BRCA1 c.5152+5G>T, identified in a HBOC family (Figure 6A). This rare variant is located in the donor splice site of intron 18 at a highly conserved position, and was predicted by several splicing in silico programs to abolish the splice site (data not shown). We obtained blood samples from patients that are heterozygous for BRCA1 c.5152+5G>T (proband and father) as well as samples from negative individuals (proband's mother and sister). To characterize the variant we performed whole transcriptome sequencing (WTS) on the proband and control blood samples. WTS detected 14 reads supporting skipping of exon 18 (Δ18, r.5075_5152del78) and 12 reads supporting the WT transcript in the proband (Figure 6B, top). In the control blood sample, only WT reads (n = 36) were detected (Figure 6B, bottom). Primers were designed in the flanking exons and RT-PCR was performed. Digital electrophoresis analysis of RT-PCR products identified a band corresponding to Δ18, exclusively in samples from the heterozygous carriers (Figure 6C). Sanger sequencing of subcloned transcripts was then performed, which confirmed Δ18 sequence (Figure 6D), indicating it results in in-frame skipping of the important BRCT functional domain of BRCA1 (19). Finally, using CloneSeq, we detected and quantified the Δ18 transcript in heterozygous individuals, which was undetectable in non-carriers (Figure 6E). CloneSeq detected 11,824 reads supporting Δ18 (r. c.5075_5152del78) and 13,268 reads supporting WT transcript in the proband (Figure 6E, left top). In the proband's father, CloneSeq detected 7,412 reads supporting Δ18, and 6,229 reads supporting WT transcript (Figure 6E, left bottom). Quantitatively, we found heterozygous individuals to express ~40% of Δ18 transcripts (Figure 6F). Of note, analysis of the LCL harboring the pathogenic alteration BRCA1 c.5152+1G>T, affecting the same donor splice site as BRCA1 c.5152+5G>T, also led to similar expression of the abnormal transcript Δ18 (Figures 6C,F). Altogether, these data were used to reclassify BRCA1 c.5152+5G>T from VUS to likely pathogenic, and therefore a clinically actionable alteration.

All authors, with the exception of SN and DT, were employees of Ambry Genetics when they were engaged with this project. The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.