Novo assay was developed to characterize SNV/INDEL, CNV and gene fusion in 483 cancer-related genes. These genes were selected based on My Cancer Genome database (https://www.mycancergenome.org), Catalogue of Somatic Mutations in Cancer (COSMIC) and other sources (18,25). Briefly, genes containing clinically important variants and genes have been reported as cancer-related were included based on a record of reimbursement in sequencing. All exons of these genes were considered which underwent hybridization-based capture from 483 cancer-related genes (Table I). For structural rearrangements detection, introns spanning recurrent fusion breakpoints were also included. Agilent's proprietary algorithm and synthetic process was used to generate the baits. The hybrid selection was done using a pool of 120-mer RNA-based baits (Agilent SureSelect) with overlap excess 3-fold for target region. All 47660 hybrid baits for catching target region constitute 2.3 Mb genomic positions, including 7011 exons and 94 introns.

Tumor specimens were collected from non-small cell lung cancer (NSCLC) and breast cancer patients at Chinese PLA General Hospital with informed consent according to the internal Review and rules of Ethics. In the very beginning of this assay, clinical samples should match several standards as follows to ensure downstream analysis. At least 10 slices of 5 μm FFPE sections or tissues with a volume of >1 was required. For each sample, hematoxylin-eosin stained slides (Fig. 1) were prepared and reviewed by a pathologist to estimate tumor purity. All samples with <50% tumor purity were marked for tumor enrichment by microdis-section to minimize contamination from normal cells (Fig. 2).

Normal cell lines harboring the population distribution of known germ line variants were mixed, and multiplexed pools with low MAF variants were used to assess and validate the limit of variant detection. First of all, to get the variants set for assessment, we sequenced 5 cell lines from the 1000 Genomes Project (26) individually and got the SNP and INDEL sites from dbSNP database (build 146) consistent with a homozygous (MAF >90%) or heterozygous (40%<MAF<60%). To estimate the INDEL detection performance, 3 additional cell lines from COSMIC database (http://cancer.sanger.ac.uk/cancergenome/projects/cosmic/) which also were sequenced individually to get the original MAF of cancer-related somatic variants in each sample. All 8 cell lines were mixed together in designed proportions, and the expected MAF of each variant was calculated on the mixed ratios (Table II). Eventually, we achieved the 2625 variants spanning a range of expected MAF (0.5–20%) and INDEL lengths (1–40 base pair, bp) as gold-standard (Table III). Cell lines obtained from Coriell Institute (http://ccr.coriell.org/) and ATCC (http://www.atcc.org/) were routinely cultured in Dulbecco's modified Eagle's media (DMEM) with 10% heat-inactivated fetal bovine serum (FBS; Invitrogen, Waltham, MA, USA) in a 75-cm² cell culture flask. The cells were seeded into cell culture flasks at a concentration of 1×10⁵ viable cells/ml and incubated at 37°C in a humidified atmosphere containing 5% CO₂.

Generally, genome DNA extracted was performed using DNeasy Blood & Tissue kit (Qiagen, Hilden, Germany). For FFPE sample special, DNA was isolated using the GeneRead DNA FFPE kit (Qiagen, Valencia, CA, USA) following the protocol. Besides the purification of high yields of DNA from FFPE tissue sections, this kit could remove deaminated cytosine to prevent false results in sequencing (27). The ratio of absorbance at 260 and 280 nm is used to assess the purity of extracted DNA, and we used the Qubit^® Quantitation Platform to quantitated DNA. A Covaris S220 focused-ultrasonicator (Covaris, Woburn, MA, USA) was used to fragment genomic DNA (500 ng) and an Agilent Bioanalyzer 2100 (Agilent Technologies) to ensure an average fragment size of 200 to 400 base pair (bp). The library preparation after fragmentation were done using instruction manual of KAPA Hyper Prep kit. The protocol included: i) repairing the DNA ends; ii) adding ‘A’ base to the DNA fragments; iii) ligating the paired-end adaptor; iv) purifying the sample using AMPure XP beads; and v) amplifying the adaptor-ligated library and purifying the sample using AMPure XP beads. Prepared library was hybridized using NCP custom designed baits as described in SureSelectQXT (Agilent Technologies) and the product was then amplified for 14 PCR cycles. The size range of the prepared library was assessed using Agilent 2100 Bioanalyzer and qualified using ABI StepOnePlus. The concentration of each library was quantified using qPCR NGS Library Quantification kit and Protocol was used to calculate the final pooling volume to sequencing. The products were sequenced using the Illumina HiSeq X platform with paired-end sequencing runs (2×150) under Illumnina recommended protocols.

Clean data were generated by data processing steps including base calling, demultiplexing and adapter trimming. All these steps were performed using Illumina HiSeq X vendor software on default parameters. We further performed our in-house software for clean data quality control (QC) which included: i) removing read pairs if any one of the two reads containing base ‘N’ >10%; ii) removing read pairs if any one of the two reads containing base with quality below Q10 >50%; iii) trimming the 3′ end of the read from the first base below Q20; and iv) removing reads shorter than 100 bp. Clean data after QC were mapped to the human reference genome (GRCh37) using BWA aligner v0.7.8 (28) with the default parameters. PCR duplicate read removal was done using Picard 1.119 (http://picard.sourceforge.net/index.html). According to the result, a sequence metric collection was generated including the number of total reads, percentage of reads mapped, on target reads number, average target coverage and percentage of target region with >200X and 1000X coverage. Before SNV and INDEL calling, local realignment was performed using Genome Analysis Toolkit (GATK version 2.7–2-g6bda569) (29,30) with default parameters and recommended ‘known sites’ in GATK best practice (https://software.broadinstitute.org/gatk/best-practices/). For SNV detection, we denote the reference allele and the coverage of each site as r and d and denote the error rate corresponding to the base calling at read i (i = 1…d) as e_i. We used a null model to explain the data in which there is no SNV at that site and all non-reference alleles to be sequencing error. The number of variant bases (k) with e_i <1e⁻³ (associated Phred-like quality score q_i>30) in each site was then given a binomial distribution. The probability under this null model was given by the following formula:

where P(X = i|d) was the probability of observing i variants in the d reads of the site. Assuming the sequencing errors were independent across reads and occurred with probability e₀ (e₀ = 1e⁻³/3) to each non-reference allele. We could obtain

The P-value was then given by P(X≥k|d) and the cut-off (P-value <1e⁻⁶) was established to eliminate random sequencing error. For INDEL detection, we simply kept variants supporting reads >10. We also employed several filters to reduce systematic errors. Empirical filters including strand bias (Fisher's exact test, P<1e⁻⁶), site median base quality (MBQ >30), site median mapping quality (MMQ >30), variant MAF (MAF >0.5%). Variants pass filters were annotated by dbSNP b146, My Cancer Genome database (https://www.mycancergenome.org) and Oncomine database v1.4.1 to get the clinical relevant information. However, cross library contamination may occur and a report would not be generated once the sample contained >10 variants with low-MAF (MAF ≤10%) in dbSNP. In the report stage, all annotated variants with MAF ≥5% would be reported and other cancer-related variants would be validated by 3dPCR. The whole workflow for the data analysis is outlined in Fig. 3. The parameters and descriptions used are listed in Table IV.

To measure the effect of our approach, we compared the pooled cell-line result with GATK, a widely used software. We followed the ‘GATK best practice’, the ‘IndelRealigner’ parameter ‘LOD_Threshold_For_Cleaning’ was 0.3, the ‘BaseRecalibrator’ was with default parameters, the SNV/INDEL calling type was ‘HaplotypeCaller’ with parameters ‘standard_min_confidence_threshold_for_emitting’ as 10 and ‘standard_min_confidence_threshold_for_calling’ as 30.

For sensitivity estimation, variants detected in pools would be assigned as true positive (TP), or false negative (FN) if not detected. Sensitivity was calculated as TP/(TP+FN). For specificity estimation, the pool variants also detected in the pure sample were assigned as true positive (TP), or false positive (FP) if none was detected. PPV was calculated as TP/(TP+FP).

dPCR is a method used in absolute quantification analysis of clonally amplified nucleic acids (including DNA, cDNA, methylated DNA or RNA). With dPCR, a sample is partitioned so that individual nucleic acid molecules within the sample are localized and concentrated within many separate regions. After PCR amplification, nucleic acids may be quantified by counting the regions that contain PCR end-product, positive reactions. Here, we used the QuantStudio™ 3D Digital PCR System platform (Life Technologies) regarding SNP mutation quantitation. For dPCR, the first step is preparing and loading samples onto QuantStudio™ 3D Digital PCR 20K chips. Mutations were analysed by TaqMan^® SNP Genotyping Assays (Life Technologies), which containing TaqMan^®-MGB probes and primers. We prepared 15 μl reaction mixes according to the manufacturer's instructions, and loaded 14.5 μl onto each chip. The Mix contains ROX^® dye, which served as a passive reference. After chips were loaded, we run the Digital PCR 20K Chips with a ProFlex™ 2x Flat PCR System under the following conditions: 96°C for 10 min, 39 cycles at 56°C for 2 min and at 98°C for 30 sec, followed by a final extension step at 56°C for 2 min. After thermo-cycling, we analyzed the prepared chips using dPCR instrument.

ARMS-PCR is a real-time PCR-based test which covers the 29 EGFR hotspots from exon 18–21. The assay was performed according to the manufacturer's protocol for the ADx EGFR29 Mutation kit (Amoy Diagnostics, Co., Ltd., Xiamen, China) with the MX3000P (Stratagene, La Jolla, CA, USA) real-time PCR system. Template DNA (0.4 μl), 3.6 μl deionized water and 16 μl other reaction components was used in the RT-PCR reaction system. PCR was performed with initial denaturation at 95°C for 10 min, followed by 40 cycles of amplification (at 95°C for 30 sec and 61°C for 1 min). The results were analyzed according to the criteria defined by the manufacturer's instructions. Positive results were defined as [Ct(sample) − Ct(control)] < Ct(cut-off).

NCP is a NGS-based clinical test for detection of somatic cancer related mutations. DNA was extracted from tumor tissues and FFPE samples, 500 ng of which was fragmented, captured using custom-designed hybridization-based biotinylated cRNA reagents and amplified via limited-cycle PCR to enrich 7,011 exons and 94 introns of 483 cancer related genes (totaling ~2.3 million sites). We used clinical samples to generate the bioinformatics pipeline for data analysis (Table IV) and cell lines to validate the whole work flow. For the 8 single cell lines, using the Illumina HiSeq X platform, achieving an average of 13,330 Mb (SD=3,995 Mb) total bases with 38.09% on-target (SD=4.78%), target regions were sequenced to 2148X (SD=537X) median coverage across targeted bases, with 99.05% (SD=0.28%) of targeted bases covered by at least 200 reads (Table V). The 2453 SNV and 172 INDEL detected in single cell line consistent with database would be used for assessment of SNV/INDEL detection. Pools of mixed cell lines were used to get the relationship between median coverage and performance, which achieved total bases of 4,762, 10,896 and 16,351 Mb, the median coverage of 1,029X, 2,237X and 3,194X (Table VI). Due to the high sensitivity NGS benefit from high coverage, the hotspot mutations with MAF <5% detected by this assay in 35 FFPE samples were confirmed by dPCR. All samples used in this test are summarized in Table VII. Finally, 33 hotspot mutations detected by NGS in FFPE samples with a MAF from 2 to 63% in NGS were tested by ARMS-PCR.

SNV detection was performed using a Binomial methodology allowing the detection of low MAF somatic mutations across the 2.3 Mb assayed with high sensitivity. For the mixed cell line pools, overall SNV detection performance was high, the results of different depth are shown in Table VIII, for an average depth of 2237, 100% (95% CI, 95.1–100%) of SNV at MAF >10% were successfully detected, as well as 99% (95% CI, 98.6–100%) of SNV at MAF 5–10%. The detection of SNV with MAF between 0.5–5% performance was 92.2% (95% CI, 90.7–93.5%) (Fig. 4A and C and Table VIIIA). In addition, high sensitivity was accompanied with good PPV (the fraction of SNV calls in the pools can also be detected in any of the individual cell lines; Table VIIIB) 99.2% (95% CI, 99–99.4%). The false positives may be due to variants with such a low MAF (<5%) no difference with sequencing noise could hardly be identified. A dPCR confirmation for cancer-related SNV with MAF <5% reported by NGS is necessary before reporting.

For INDEL detection, we simply discarded the variants supporting less than 10 reads. The results of different depth are shown in Table IX, for an average depth of 2237, 100% (95% CI, 29.2–100%) of INDEL at MAF >10% were successfully detected, as well as 94.7% of INDEL (95% CI, 74–99.9%) with MAF between 5–10%. Low MAF sites detected performance was 91.5% (95% CI, 85–100%), the performance of variants with MAF <0.5% was also calculated (Fig. 4B and D and Table IXA). Few false-positive calls were observed, with a PPV of 98.2% (95% CI, 97.2–98.9%) (Table IXB). Like SNV detection, due to the false positive under 10%, a dPCR confirmation of these cancer-related INDEL with MAF <10% before reporting is needed.

We evaluated the performance of our bioinformatics pipeline with the cell line models above, focusing on two key steps of our approach. First, we applied statistical models that allow for the identification of a mutation at low MAF from random errors in Illumina sequencing. Second, we used priori knowledge to identify systematic errors always accompanied with specific characteristics, such as strand bias and low base/mapping quality. To measure the effect of our approach, we compared the pooled cell-line result with GATK - widely used software. The GATK detection sensitivity of SNV with MAF >10% was 64.38% (95% CI, 52.3–75.3%), and SNV with 5%<MAF<10% was under 10% but the PPV was 100% (95% CI, 99.7–100%). The sensitivity of INDEL with MAF >10% was 67% (95% CI, 9.4–99.2%), and a high PPV 100% (95% CI, 99–100%) (Tables X and XI), possibly because this widely used tool is designed for whole-genome or whole-exon sequencing data with relatively low depth and variants with high allele frequency, which underline that appropriate filters for ultra-deep sequencing data analysis were critical. Actually, compared with slight performance upgrades under increased coverage depth, the effect of appropriate filters was remarkable in this test.

The above studies demonstrate that the NGS-based test has the performance characteristics necessary to accurately detect SNV and INDEL. We further validated test accuracy by comparisons to dPCR for 35 FFPE cancer specimens. To assess the accuracy of low MAF SNV and INDEL detection in routine clinical cancer samples, we selected 35 FFPE resection specimens (31 non-small cell lung cancer, 1 parathyroid carcinoma, 3 breast cancers) previously tested for hotspot mutations in PIK3Ca, EGFR, KRAS and BRAF by NGS, every hotspot mutations detected by NGS, but with MAF <5% would be tested by dPCR. In addition, 32 of 35 (PPV=91.43%, 95% CI, 76.94–98.20%) variants have been supported to be true-positive by dPCR (Tables XII and XIII). Three variants were present at <3% MAF in NGS that were not detected by dPCR. The detected MAF of the two technologies is shown in Fig. 5. Finally, we random selected 33 FFPE samples (NSCLC) with hotspot mutations and performed the ARMS-PCR to verify the overall PPV of our assay. As a result, all 33 mutations could be detected by ARMS-PCR and the PPV was 100% (95% CI, 89.42–100%; Table XIV).