Performance of amplicon-based next generation DNA sequencing for diagnostic gene mutation profiling in oncopathology

Study population

For this study, a series of FFPE tumor tissue–derived DNAs from the files of the Department of Pathology (VU University Medical Center (VUmc), Amsterdam, The Netherlands) was used. First, a learning (n = 58) and test (n = 45) FFPE-DNA series was selected, with 78 % (45/58; learning set) and 91 % (41/45; test set) of tumor samples harboring known somatic mutation(s) in EGFR, KRAS, BRAF and/or PIK3CA as identified during routine diagnostic practice (period March 2009 to August 2013). Subsequently, a prospective series of FFPE-DNA from 75 tumor samples that were routinely screened for somatic mutations in EGFR, KRAS, BRAF and PIK3CA was used (period September 2013 to December 2013). In all specimens, the mutational status of EGFR exon 18–21, KRAS exon 2/3, BRAF exon 15 and PIK3CA exon 9/20 was assessed by HRM assay followed by Sanger sequencing of HRM-PCR products with an aberrant melt curve, essentially as described previously [25–27]. Primers used for PIK3CA exon 9 were 5’-GTAAAACGACGGCCAGGACAAAGAACAGCTCAAAGCAA-3’ (forward primer) and 5’-CATTTTAGCACTTACCTGTGACTCCA-3’ (reverse primer) and for PIK3CA exon 20 were 5’-GTAAAACGACGGCCAGTGGAATGCCAGAACTACAATCTT-3’ (forward primer) and 5’-TGTGTGGAAGATCCAATCCA-3’ (reverse primer), respectively. The analytical sensitivity of HRM-sequencing is ≤5 % of mutant DNA in a wild-type background [25–27]. All samples were used in compliance with the respective institutional ethical regulations for surplus material [28]. Specimens were processed according to a routinely applied protocol [25–27] in which the FFPE tumor tissue was manually macro-dissected from serial paraffin sections guided by a H&E stained tissue section on which the tumor was marked by a pathologist. DNA was subsequently isolated using a QIAamp DNA FFPE Tissue kit according to the instructions of the manufacturer (Qiagen, Venlo, The Netherlands). DNA concentrations were measured using a Qubit 2.0 Fluorometer (Life Technologies (Carlsbad, CA). DNA isolated from human cancer cell line HCT116 (KRAS p.G13D; PIK3CA p.H1047R; CTNNB1 p.delS45) was used as a mutation control, and an XY reference DNA pool (Kreatech Diagnostics, Amsterdam, The Netherlands) was used as a normal, wild-type, control.

TSACP-MiSeq-NGS

The TruSeq Amplicon Cancer Panel (TSACP; Illumina Inc, San Diego, CA, USA) comprises 212 amplicons from 48 genes that are simultaneously amplified in a single-tube reaction. We used 5 μl of each FFPE-DNA sample (total of 150 ng DNA, unless otherwise specified) as input for amplicon library preparation according to the manufacturer’s instructions. The resulting sequence library was quantified using a Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA). All sample libraries, including those of samples with low yield, were used for subsequent sequencing. Up to 24 differently barcoded, individual sequence libraries were equimolarly pooled prior to sequencing. This multiple-sample sequence library pool was loaded on a MiSeq Personal Sequencer (Illumina) with a molarity of 10 pM, and a PE150 run was performed using a MiSeq Reagent Kit v2 (300 cycles) (Illumina) according to the manufacturer’s instructions. TSACP-MiSeq-NGS generated consistent variant allele frequencies (VAF) in a 10-fold testing of HCT116 DNA (median VAF KRAS p.G13D of 51 % (range 48–52 %); median VAF PIK3CA p.H1047R of 49 % (range 48–51 %); median VAF CTNNB1 p.delS45 of 50 % (range 48–52). TSACP-MiSeq-NGS detected KRAS p.G13D, PIK3CA p.H1047R, and CTNNB1 p.delS45 in HCT116 cell line serial dilutions down to ≤5 % VAF [29].

Bioinformatics pipeline

First, the sequence reads were de-multiplexed by Bcl2fastq (Illumina), cleaned for adapter sequences by Cutadapt [30] and trimmed for trailing nucleotides with a quality lower than Q30 by the FASTX-Toolkit [31]. Paired-end reads were merged into a single read using overlapping sequences with Flash [32]. Next, the generated single reads were aligned to the human reference sequence build GRCh 37/hg19 by BWA [33]. The assessment of sequence variants and indels in the single reads, relative to the human reference sequence, was performed using an in-house developed variant caller tailored to FFPE-DNA. The variant caller utilizes an error model which takes into account errors related to the sequencing technology (E ^sbs) and formalin fixation (E ^fix), as the latter may introduce a bias in the type of nucleotide substitutions observed [11, 21, 23]. These errors are estimated by

where N ^mismatch is the number of nucleotides in the reads which do not match the reference genome and N ^match is the number of nucleotides matching the reference genome, and

where N ^i>j is the count of i to j nucleotide transitions and N ⁱ is the occurrence of nucleotide i in the reference genome. The error model is used with the binomial distribution to calculate a PHRED scaled quality score (QUAL) for each variant as

where Nvar ¹ is the number of variants at locus l , and D ¹is the depth at locus l.

Annotation for dbSNP entries and mutation effects were performed by SnpEff [34]. Finally, by using web standards (HTML) a report was generated that summarizes per specimen the results for all amplicons, including the sequence depth and, if present, the called mutation (both on DNA and protein level), VAF, QUAL and dbSNP identifier. The Integrated Genomics Viewer [35] was used for visual evaluation of the data. The bioinformatics pipeline is available at https://bitbucket.org/daoud/falco.

Statistical analysis

TSACP-MiSeq-NGS was performed blinded for the HRM-sequencing data, and the results were compared afterwards. Results of both assays comprised singular measurements. Variants that were designated as common and without known clinical association by dbSNP and/or the 1000 genomes project [36] were excluded from analysis. The level of agreement between the assay findings of the two different methods was depicted as the percentage of concordant findings. In this analysis, only nucleotide substitutions that could be detected by both assays in exons of the EGFR, KRAS, BRAF and PIK3CA genes were considered. Results were considered concordant in case of complete genotype agreement between both assays, and discordant in case of genotype disagreement. The level of agreement between both assays was determined using kappa statistics. Specimens for which TSACP-MiSeq-NGS failed to generate a result were excluded from downstream statistical analysis. Analyses were performed using SPSS statistical package version 20 (SPSS Inc., Chicago, IL, USA). The level of statistical significance was set at 0.05.

Performance of TSACP-MiSeq NGS in a retrospective tumor learning/test set

We first evaluated a series of 58 FFPE-DNA samples of routine diagnostic tumor specimens with known somatic mutations in EGFR (exon 19 deletion, n = 10; exon 21 point mutation, n = 6; exon 20 insertion, n = 4; exon 20 point mutation, n = 2; p.T790M, n = 5), KRAS (codon 12 point mutation, n = 14; codon 13 point mutation, n = 4; codon 61 point mutation, n = 4), BRAF (codon 600 point mutation, n = 2; codon 597 point mutation, n = 1) or PIK3CA (codon 1047 point mutation, n = 1) as determined by HRM-sequencing. These mutation data were used as a reference to tailor the NGS methodology and bioinformatics pipeline for reliable use of FFPE-DNA. Three critical minimum parameter settings were deduced from the learning set to consider an FFPE-DNA sample qualified for TSACP-MiSeq-NGS, namely (i) a minimal sequence library yield of 1 nM, (ii) a minimal amount of 100.000 reads, and (iii) a minimum of 90 % of amplicons with minimally 100 times coverage. Samples that did not meet these quality criteria (n = 6) exhibited a marked underrepresentation of mutations in the NGS data (i.e., 5 out of 7 mutations not identified), supporting that these FFPE-DNA specimens should be disqualified for diagnostic mutation analysis by TSACP-MiSeq-NGS. Critical minimum settings for scoring a mutation among qualified specimens were (i) a minimal VAF of 10 %, (ii) a minimal amount of 100 reads at the given position, and (iii) a minimum QUAL score of 990. Using these thresholds, 50 out of 52 (96 %) qualified FFPE-DNA specimens showed similar assay findings for KRAS, PIK3CA, BRAF and EGFR as HRM-sequencing. When lowering the VAF threshold to 5 %, the remaining two qualified FFPE-DNA samples also yielded mutations by TSACP-MiSeq-NGS (i.e., point mutations in KRAS and PIK3CA, both at 9 % VAF), which were concordant with the HRM-sequencing outcome. However, the use of this threshold for TSACP-MiSeq-NGS was at the cost of calling various additional alterations in the learning set that could not be confirmed by HRM-sequencing, although they were within the detection range of the latter technique, including BRAF p.A598G and EGFR p.V765L. As sequence artifacts may mimic clinically important mutations, we considered the latter unfavorable for diagnostic mutation analysis and, therefore, chose to fix the limit of detection of TSACP-MiSeq-NGS to 10 % VAF to ensure reliable application with FFPE-DNA, at the expense of sensitivity.

The trained minimum parameter settings with chosen thresholds for quality and mutation calling were subsequently validated using FFPE-DNA samples from an independent series of 45 tumors representing various somatic mutations in KRAS, PIK3CA, BRAF or EGFR as determined by HRM-sequencing, including single or double nucleotide substitutions and indels. The sample details and sequencing statistics are listed in Table 1. Eight of the FFPE-DNA samples (18 %) did not pass the quality criteria for diagnostic mutation analysis by TSACP-MiSeq-NGS. These included all samples with a DNA input of 30 ng or less (n = 4), showing that these low DNA amounts were insufficient for TSACP-MiSeq-NGS. In the remaining 37 qualified specimens, a full concordance in assay findings of TSACP-MiSeq-NGS and HRM-sequencing for KRAS, EGFR, BRAF and PIK3CA was observed in 95 % (35/37; kappa = 0.92; 95 % CI = 0.81–1.03). In two samples (V23, V29), the KRAS mutations detected by HRM-sequencing were not called by the TSACP-MiSeq-NGS approach. Lowering the thresholds would have identified the respective mutations (V23, amplicon KRAS2, 1475 reads, p.Q61L, 8 % VAF, QUAL 1000; V29 amplicon KRAS7, 2540 reads, p.G12D, 6 % VAF, QUAL 840), but again at the cost of various alterations that could not be confirmed by HRM-sequencing. One sample (V6) with somatic mutations in codons 790, 833 and 835 of EGFR as observed by HRM-sequencing, was only called for codon 790 as the other codons of this gene are not covered by the amplicons in the TSACP.

Performance of TSACP-MiSeq NGS in a prospective tumor series

FFPE-DNA samples from 75 consecutive tumors for which mutation analysis was requested in a routine pathology setting, were evaluated by TSACP-MiSeq-NGS according to the validated test algorithm. Of these, 5 (~7 %) yielded limited DNA (maximum input <30 ng). These were not further considered for TSACP-MiSeq-NGS. Seven (~10 %) of the remaining specimens did not meet the quality criteria for diagnostic mutation analysis by TSACP-MiSeq-NGS (i.e., 4 showed an insufficient library yield and 3 an insufficient read count and/or pass rate). Altogether, this translates into a success rate of TSACP-MiSeq-NGS of 84 % (63/75) in the routine pathology setting used. Among the 63 successfully sequenced specimens, we obtained an average number of 781,871 reads (range 122,578– 3,227,821), and an average percentage of passed amplicons (>100 × coverage) of 98 % (range 90–100 %), with failed amplicons mainly from the CDNK2A gene. The genes with mutations (not including common SNPs) called in these specimens are depicted in Fig. 1. A full concordance for KRAS, EGFR, BRAF and PIK3CA mutations that could be detected by both assays was observed in 95 % (60/63) of the specimens (kappa = 0.92; 95 % CI = 0.84–1.01). A total of 4 mutations in KRAS and EGFR found in three specimens by HRM-sequencing, respectively, were not called by TSACP-MiSeq-NGS (i.e., case D33 KRAS p.G12V, case D58 EGFR p.L858R and p.T790M, and case D59 KRAS p.G12D). TSACP-MiSeq-NGS did not identify alterations other than those detected by HRM-sequencing in these target regions. Nonetheless, extra mutations in other cancer-related genes, such as TP53 in 35 % (22/63) of cases, were clearly identified (Fig. 1).