Validation of amplicon-based next generation sequencing panel for second-tier test in newborn screening for inborn errors of metabolism

Validation samples

The validation samples consisted of two reference genomes, RM 8398 (ethnic Utah/Mormon) and RM 8393 (ethnic Chinese) (purchased from NIST, Gaithersburg, USA). A total of three validation experiments were performed, each of which included a 16-plex sequencing run, spiked-in with 2% PhiX DNA as technical control. One sample was triplicated in each run for validation of technical performance on repeatability and between run repeatability. A total of 3 × 16 = 48 sequencing results were obtained. The reference genome samples were used for assessment of sensitivity and specificity.

Panel design

A custom AmpliSeq panel (Illumina, San Diego, USA) was designed for this study, wherein was amplicon-based, targeting the exons and ±20 base pair (bp) of intron-exon-boundaries of 87 genes and 289,264 bps in the human genome with 2,263 amplicons (Panel design summary shown in Table 1). The primers were separated into two pools for library amplification. The panel could detect genetic variants for 87-IEM diseases that are potential targets of NBS or for other clinical purpose (Table 2). For reference, a Browser Extensible Data (BED) file containing the genomic positions of the missed bases in exons is provided as Supplementary Data 1, and the exons containing missed bases in each gene were summarized in Supplementary Data 2.

DNA extraction for DBS

For DBS collected on Whatman 903 Proteinsaver Cards (Cytiva, Vancouver, Canada), six 3.2 mm diameter DBS discs should be used for DNA extraction. Extraction was performed manually using column-based QIAamp DNA Mini Kit (QIAGEN, Germantown, USA) with modification using manufacturer-supplied buffers ATL, AL, AW1, AW2, EB, and proteinase K. Three-hundred and sixty microliters buffer ATL for tissue lysis was added to the punched DBS discs, followed by incubation at 85 °C for 10 min; 40 μL proteinase K stock solution was added and incubated at 56 °C for 1 h; 400 μL buffer AL was added and incubated at 70 °C for 10 min. Four-hundred microliters ethanol (96–100%) was added to the mixture before applying to the QIAamp Mini spin column and centrifuged at 6,000×g for 1 min with filtrate discarded, followed by washing with AW1 and AW2. DNA was eluted from column with the same 50 μL buffer EB (QIAGEN, Germantown, USA) twice with 1 min incubation at room temperature, followed by centrifugation at 6,000×g for 1 min. The extracted DNA quantity and quality were accessed by Qubit DNA High Sensitivity (HS) Assay Kit with Qubit 2.0 Fluorometer and NanoDrop spectrophotometer (ThermoFisher Scientific, Waltham, USA) respectively.

Library preparation and targeted sequencing

Library preparation was carried out using AmpliSeq Plus Library Prep Kit (Illumina) following the manufacturer’s protocol. About 7.5 ng of DNA was used to synthesize library. AmpliSeq UD indexes were used for indexing. At the end, 27 μL of eluted library DNA was collected to a DNA LoBind tube (Eppendorf, Hamburg, Germany). NGS libraries were quantitated by Qubit dsDNA HS assay. Library size distribution was checked by Agilent High Sensitivity DNA Kit on 2,100 Bioanalyzer (Agilent Technologies, Santa Clara, USA), from which a distinct peak between 300 and 400 bp should be observed, with negligible level of contaminating DNA species. Each sequencing library was diluted to a concentration of 1 nM with buffer EB and 16 libraries were pooled together. The pooled library was then diluted to 90 pM and spiked with 2% 90 pM PhiX DNA as technical control. The pooled library was then loaded to the sequencing cartridge and subject to massive parallel sequencing to generate 150 bp pair-end reads in Illumina iSeq-100 System.

Data analysis

The iSeq-100 built-in software, DNA Amplicon, generated FASTQ output files from the raw reads, performed initial alignment and determined QC metrics. Sequence Analysis Viewer (SAV; Illumina) and FastQC (Bioinformatics Group, Babraham Institute, Cambridge, UK) apps were also used to monitor sequencing QC metrics of pooled and individual library, respectively.

We employed a commercial pipeline, the NextGENe software (v2.4.2.3; Softgenetics, State College, USA), to perform read alignment against the reference genome (GRCh37/hg19), filtering and variant calling. Data analyses targeted all the coding exons ±20 bp to include the splice site intronic region. For variant calling, the criteria included >20% mutation frequency, >3 allele count and >5 total coverage count; homozygous variants were exempted from these cutoffs; in-read phasing (merging) of adjacent variants with a maximum gap of 1 bp and phaseable read percentage of >50% was allowed. Additional custom filtering criteria was imposed to minimize false-positive rates by using a proprietary quality score (cutoff = 4) which was composed of adjusted coverage score, mismatch score and wrong allele score. A variant was regarded as valid only if a minimum of 20 reads at its position was observed. A mutation report (VCF file) was created annotating all variants. Human Gene Mutation Database (HGMD) Professional (QIAGEN, Germantown, USA), Alamut Visual (Sophia Genetics SA, Saint-Sulpice, Switzerland) and Integrative Genomics Viewer (IGV) [13] were used for further characterizing the variants and raw reads. Single nucleotide variant (SNV) and small indel were targeted for curation.

External quality assurance (EQA) program

We joined the European Molecular Genetics Quality Network (EMQN) pilot EQA scheme for NGS (Germline; 2020) as EQA. The provided DNA sample was sequenced in another 16-plex run, and the BED file of the original panel design covering the whole-panel, instead of coding exons only, was given to the EQA program to analyze the data.

DNA extraction, library preparation, and sequencing metrics

All libraries showed satisfactory yield and quality, i.e., concentration ranged from 2.67 to 11.6 ng/μL and average size ranged from 359 to 390 bp with a distinct and dominant peak between 300 and 400 bp. The total yield of each run ranged from 1.74 to 1.92 giga bps. The overall sequencing quality of the three runs was good with %Q30 between 91 and 92% (Table 3). PhiX alignment rate is between 1.3 and 1.4% (compared with the expected 2%; Table 3). All individual library passed the FastQC metrics except one library showed distinctively high duplication level. This sample also showed a low uniformity (defined as % of base with >0.2× mean coverage) thus regarded as fail.

Coverage and repeatability

After alignment against the reference genome (GRCh37/hg19), on average one library got about 640,000 reads on targets, with an average coverage of 208× ±20 (1 SD) (Table 3). For intra-run repeatability, the triplicated sample had mean coverage of 195, 206, and 223 in each run, corresponding to CV of 9.1, 10.5, and 5.2% respectively. Combining the above to calculate inter-run repeatability, the average coverage is 208× with a CV of 10.1%.

For each of the 87 genes, we looked at the percentage of >20× coverage in reportable region from the three runs (Table 2). Seventy-five genes achieved an average of ≥90% >20× coverage, of which 8 genes achieved 100% >20× coverage in all three runs. Another eight genes had between 80 and 90% >20× coverage and were also useable for reporting. CYP21A2 was disqualified due to its very high homology with its pseudogene CYP21A1P. Another gene, IDUA, has low >20× coverage at 55.9% and will not be used for reporting. NAGS and GCSH had the next lowest >20× coverage at around 70% but were still considered useable for reporting.

Sensitivity

We used the variants in the reference genomes RM8398 (ethnic Utah/Mormon) and RM8393 (ethnic Chinese) to evaluate the sensitivity of the panel by looking at the proportion of true-positive variant called. The sensitivities were 90.9% (442 out of 486 variants; SNV and indel combined) and 91.4% (367 out 401 variants; SNV and indel combined) respectively which were considered acceptable. Specificity was calculated based on the proportion of false-positive over the total number of base in the reporting region achieving >20× coverage, so that specificities were 99.94% (98 false-positives out of 163,932 bases) and 99.95% (87 false-positives out of 163,932 bases) for RM8398 and RM8393 respectively. Many of the false-positives were found near the 3′ end of the mapped reads with lower coverage and quality score, probably due to a drop in sequencing quality near the 3′ end. These assay specific artefacts were recurrently observed in all samples, and could be easily identified and avoided using genome browser during NGS data reporting.

EQA

The results returned from EMQN pilot EQA scheme for NGS (2020) was satisfactory (Table 4), with the median depth at called variants being 178×. The sensitivity for SNV and indel were 97.89 and 81.58% respectively. The precisions were relatively lower for both SNV and indel at 77.44 and 53.4% respectively, due to the fact that the whole-panel BED file was submitted for analysis, which includes more non-coding regions and 3′ end of reads. The latter are known to more likely exhibit lower sequencing quality and PCR artefacts. The artefacts at 3′ end are also more likely to occur in tandem hence making indel more likely to occur than SNP. Fortunately, these artefacts can be easily identified using genome browser. Since we were the only participant who used iSeq and the custom AmpliSeq panel (out of 293 participating labs; 77.6% of which used one of the Illumina NGS platforms), our results cannot be directly compared with the peer group. Our EQA results show good sensitivity and acceptable precision, both of which can be improved if filtering criteria is optimized and only coding exons ±20 bp are analyzed.

Turnaround time (TAT)

For a 16-plex run, DNA preparation from DBS can be finished in half a day (Day 1). Library synthesis and pooling can be accomplished within 1.5 days (Day 1 and 2) followed by overnight sequencing for about 22 h. Bioinformatic analysis and QC procedure can be done within 3 h at Day 3. With normally only a few variants per gene being discovered and one to two genes per sample needing to be interpreted, an optimized workflow should allow reports to be ready in about 2.5 working days.

Optimizing sample preparation and scrutiny

DNA extraction was performed manually using column-based method; in our experience, the minimum yield was about 60 ng, which is more than required for quantitation and library construction (about 12 ng). The excess DNA should be enough for repeat or other tests such as independent single nucleotide polymorphism (SNP) genotyping assays. Considering the urgency of NBS reporting, automated DNA extraction with magnetic particles could be an alternative method with the advantage of better speed, yield, purity and consistency, as well as reduced chance of sample swapping. The complicated sample and library preparation workflow in targeted NGS has introduced the possibility of sample-mix up or cross contamination [14, 15]. Therefore, procedures to scrutinize sample identity are important when implementing the test in order to minimize the chance of misreporting. Possible way includes in silico identity check to look up a set of highly heterogeneous SNP genotypes and sex prediction from the sequencing reads; the amalgamated SNP genotype should not be identical among a batch of 16 samples and the predicted sex should match with the reported sex. Furthermore, independent PCR-based SNP genotyping assays can be performed and the result compared with that from the in silico SNP identification pipeline, in order to minimize the chance of sample-mixing/swapping.

Second-tier NBS test through targeted NGS

The introduction of NGS would be beneficial to current NBS practice. First-tier MS/MS and biochemical approaches allow detection of multiple IEM including amino acid, organic acid and fatty acid disorders. However, these tests may have short-comings of relatively high false-positive or false-negative rates. This might be due to factors such as prematurity, total parenteral nutrition, maternal defects and the age of disease onset, which affect metabolite concentration and lead to false results. Though two-tier testing might help to screen out suspected cases for specific IEM, limited applications are available. To complement with the current NBS approach, targeted NGS could be implemented as second-tier test. The drawbacks of first-tier biochemical tests could be addressed through the integration of targeted NGS, to minimize the false-positive or negative results due to borderline metabolites concentration.

According to the statistics of Hong Kong NBS since 2015, several conditions are found to have relatively high incidence of false-positive or false-negative results [4], therefore, these conditions would be excellent candidates for second-tier test using our assay. For example, our NGS assay could be used as second-tier test for samples with borderline range of free carnitine level (e.g., 5.0–6.5 μmol/L). The absence of pathogenic variants in the SLC22A5 gene, the causative gene of systemic primary carnitine deficiency [16], could reduce false-positive cases.

On the other hand, false-negative cases were frequently reported for citrin deficiency in Hong Kong [4] and elsewhere [17, 18]. By adopting second-tier genetic testing, we can extend the detection of patients with citrin deficiency with a relatively mild citrulline elevations, higher than the upper reference limit but below the NBS reporting cutoff (e.g., citrulline level of 25–35 μmol/L).

Limitations

Short-read mapping is not suitable for genes with highly homologous genomic regions such as pseudogenes [19], in this case CYP21A2 is interfered by its pseudogene CYP21A1P. While the assay can pick up most of the pathogenic variants in the exonic and splice region, other variants within the promoter regions, deep intronic regions, or regulatory elements outside of the targeted regions cannot be detected. Also, gross insertions and deletions within the target regions could not be detected by this assay at the time of reporting. These variants could possibly be detected using other optimized post-analytical pipeline specifically designed for copy number variant (CNV) detection. Amplicon-based panel might also be susceptible to allele drop out due to the presence of SNVs that cause defective primer annealing in primer binding sites [20, 21] or secondary structure formation in amplicon non-primer site [22], culminating in insufficient coverage and false-negative result.