Introduction

Over the last decade, the field of molecular genetic diagnostic has undergone tremendous changes. Introduction of next generation sequencing (NGS) enabled replacement of single gene tests with comprehensive gene panels, whole exome sequencing (WES) and whole genome sequencing (WGS) to increase the likelihood of identifying causal variants while decreasing the number of tests [1–4].

Currently, there is a growing clinical use of WES and WGS to identify causal variants and to discover new disease related genes in patients. However, routine sequencing of large numbers of WES or WGS remains expensive for clinical use. Several issues such as incidental findings that pose significant ethical problems, cost effectiveness and technical challenges of clinical interpretation of enormous numbers of genetic variants remain challenging [5, 6]. Moreover, WGS generates a huge amount of data that require complex bioinformatic analysis tools for data handling and storage [1].

To fully leverage the power of NGS in a large number of samples in a cost- and time-effective manner, several targeted enrichment approaches are available [7, 8]. Many diagnostic laboratories are implementing targeted enrichment NGS methods to focus on specific genes panels, or genomic regions for genetically heterogeneous diseases [9–12]. The most commonly used custom-enrichment approaches are based on capture by hybridization, PCR-based methods (highly multiplexed PCR) and capture by circularization [7, 13–17].

Neuropathic pain is a common feature of peripheral neuropathy that imposes a significant impact on patients’ quality of life and health care costs. Millions of individuals (7–10% of the general population) worldwide suffer from neuropathic pain [18]. However, not all individuals with peripheral neuropathy develop pain, and it is not possible to predict who is more or less susceptible among those with similar risk exposure [19]. Current inability to identify high-risk individuals hinders development and application of therapies to counteract neuropathic pain and to address targeted prevention strategies. Pathogenic variants in voltage-gated sodium channel (VGSC) genes expressed in the peripheral nociceptive pathway such as SCN9A, SCN10A, and SCN11A have been reported to play a key role in neuropathic pain [20–23].

In order to provide a genetic diagnosis in patients with neuropathic pain, identification of causal variants in genes encoding sodium channel subunits as well as other genes involved in painful neuropathy pathway is needed. Resolving the genetic architecture of painful neuropathy will lead to better disease management strategies, risk stratification, and counselling. Therefore, the aim of this study was to develop, validate and implement a reliable technique to rapidly and accurately re-sequence multiple genes in a large cohort of neuropathic pain patients at low cost. We present here the assessment of Molecular Inversion Probes-Next generation sequencing (MIPs-NGS) and TruSeq^® Custom Amplicon-Next generation sequencing (TSCA-NGS, Illumina, Inc., San Diego, CA, USA) methods, using a custom gene panel, to identify genetic variants in patients with neuropathic pain. For both methods, we constructed a targeted enrichment kit to capture the coding and exon-flanking intron sequences of nine sodium channel genes (SCN3A, SCN8A-SCN11A, and SCN1B-SCN4B). We applied the two methods to test 166 different patients and systematically compared the sensitivity, specificity, flexibility, targeting efficiency, reproducibility of performance and cost effectiveness of MIPs-NGS and TSCA-NGS approaches.

Materials and methods

Targeted enrichment, paired-end sequencing and data processing

Targeted enrichment kits were constructed for this study to capture the coding and exon-flanking intron sequences (±20 base pairs [bp]) of nine VGSC genes; SCN3A, SCN8A-SCN11A and SCN1B-SCN4B (data in S1 Table). The probes were designed for the two methods using their respective informatics pipelines and methods are provided in the subsections below. Probe features and sequencing characteristics are given in Table 2.

MIPs-NGS.

Two hundred and seventy-six molecular inversion probes (MIPs) were designed to capture 37,467 bp that represent all exons and exon-flanking intron sequences of nine VGSC genes. A modified version of MIPgen tool (http://shendurelab.github.io/MIPGEN/) was used to design MIPs. Data in S1 Fig shows a flowchart of the MIPs design using MIPgen, including required MIPgen criteria.

All probes were fixed in to 77-80-mer in length and each MIP contains two targeting arms; an extension arm ranging in length from 16 to 20 nucleotides (nt) and a ligation arm ranging in length from 20 to 24 nt. These arms were joined by a 30 nt common linker sequence which contains two universal PCR primer sites (complimentary probe arms are available upon request). MIPs were designed to have an overlapping of 20 bps and every base in the targeted region should be covered and captured at least by one probe. All MIPs with high arm copy count (>5x) were excluded. Furthermore, 3’ and 5’untranslated regions (3’ and 5’ UTRs) were not included in the design. Common single nucleotide polymorphisms (SNPs) (>1%) in the extension and ligation arms of the MIPs were excluded whenever possible. In some cases where SNPs were not avoidable, two MIPs were designed to match both wildtype and variant genotype for the same locus (n = 12). Constructed BED file of the targeted regions and MIPs corresponding to these nine candidate genes were uploaded to UCSC genome browser (http://genome.ucsc.edu/) in order to validate each designed MIP.

To reduce costs, the standard gap-fill length of 112 nt between the extension and ligation arm (region of interest) of the MIPs [24] was adapted to 220–230 nt. Each MIP contains a unique molecular identifier (UMI) of 5 nt to remove duplicates introduced by PCR amplification and sequencing. Probes were synthesized by Integrated DNA Technologies (IDT, Iowa, USA) and delivered individually. Subsequently, probes were equimolarly pooled and phosphorylated at 5’ end of the probe (data in S1 File).

Experimental workflow was done by following standard protocols [16, 24, 25, and data in S1 File and S2 Fig]. In brief (Fig 1A), 50–100 ng of high quality, non-fragmented genomic DNA was used for hybridization. After gap filling and ligation, circularized DNA molecules were used as template in PCR with universal primers complementary to the linker sequence. Then, sample-specific barcode sequences and Illumina adaptors were introduced during the PCR amplification step. Fragment size, quality, and quantity of the amplified captured material were determined by Qubit^® 2.0 Fluorometer using the Qubit^® dsDNA BR assay kit (Life technologies, Bleiswijk, The Netherlands) and Agilent 2100 bioanalyzer using the Agilent DNA 1000 Kit (Agilent Technologies, Santa Clara, CA, USA) following manufacturer’s instructions. Next, samples were pooled and purified using Ampure XP beads (Beckman Coulter, Inc, Brea, California) according to manufacturer’s instructions. Pooled samples were then paired-end sequenced (2 x 150 bp or 2 x 250 bp) using the MiSeq or NextSeq 500 Instrument (Illumina, Inc., San Diego, CA, USA) to achieve >30x coverage per bp.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. MIPs-NGS and TSCA-NGS workflows summary.

Graphical depiction of MIPs-NGS (A) and TSCA-NGS (B) methods for library construction and capture protocol for targeted next-generation sequencing.

https://doi.org/10.1371/journal.pone.0238467.g001

After a pilot experiment with 10 DNA samples, MIPs with poor performance (15–30 reads) were rebalanced by adding a 5 fold concentration of each poorly performing probe to the original MIPs pool, and MIPs with no sequence reads were replaced by adding new MIPs to the original MIPs pool. Then, the pilot experiment was repeated as described above.

Sequenced data was analyzed by using our in-house MIPs-targeted NGS data analysis pipeline. This pipeline aligns sequenced data to the human reference sequence GRCh37, trims probe arm sequences, de-duplicates the data on the basis of UMIs, and annotates variants according with information and frequencies from ExAC, dbSNP, cadd, Gencode, and calculates coverages per sample, number of bases, mean and median coverage per MIP target [26].

Due to the coverage depth used for this approach, no copy number variant (CNV) analysis has been performed.

TSCA-NGS.

One hundred and eighty-five oligonucleotides were designed to capture 38,258 bp that represent all exons and exon-flanking intron sequences of nine VGSC genes. The design of oligonucleotides was assessed by the free tool Design Studio, which is available on Illumina website (https://designstudio.illumina.com/ (Illumina Inc., San Diego, CA, USA)). Targeted regions were selected by entering the genes name to the tool and only coding exons were selected. All probes were fixed in to 50-mer in length (including the adapters) and the gap-fill length was set to 425 bp (amplicons size). In case that the tool failed to locate automatically the probes hybridization sites, due to the high CG content or high SNPs rates, we manually added the genomic regions by entering the chromosome coordinates. The a priori calculation of the target regions coverage using this tool provided better results when adapting 425 bp amplicon size, instead of the 250 bp option. Increasing the amplicon size resulted in costs reduction in terms of number of oligonucleotides probes spanning the entire target region. One pair of oligonucleotides was designed for each amplicon and provided as a pool in a single tube.

TruSeq Custom Amplicon libraries preparation was performed according to manufacturer’s protocol (TruSeq Custom Amplicon Library Preparation Guide, Part # 15027983 Rev. C). In brief (Fig 1B), 250 ng genomic DNA was hybridized to the custom-made pool of oligonucleotides specific to the targeted regions of interest. Then, unbound oligonucleotides were removed by washing the samples on the provided filtering plate. Subsequently, DNA polymerase was added to extend from the upstream oligonucleotide through the targeted region, followed by ligation to the 5’ end of the downstream oligonucleotide using DNA ligase. The extension-ligation resulted in the formation of products containing the targeted regions of interest flanked by sequences required for PCR amplification. During amplification, the extension-ligation products were amplified using indexed-primers that added sample multiplexing index sequences (i5 and i7) as well as common adapters required for cluster generation (P5 and P7). After purification with PCR Clean-Up AMPure XP beads (Beckman Coulter, Inc.), the library was normalized with Library Normalization Beads. Finally, equal volumes of normalized library were combined, diluted in hybridization buffer, and heat-denatured for sequencing on the MiSeq instrument (paired-end sequenced 2 x 250 bp) (Illumina, Inc., San Diego, CA, USA). No probe dosage adjustment was required for TSCA-NGS.

To identify sequence variations in SCN3A, SCN8A-SCN11A, SCN1B-SCN4B by TSCA-NGS, Raw Fastq files were aligned to reference sequence GRCh37 using CLC Genomics Workbench (CLCbio, Qiagen, Hilden, Germany). Variant calling, annotation and coverage analysis was performed using CLC Genomics Workbench (CLCbio, Qiagen, Hilden, Germany). No CNV analysis was carried out for this approach.

Cross-comparative analysis of genetic variations detection between MIPs-NGS and TSCA-NGS

A cross-comparative analysis was performed in order to compare the sensitivity, specificity and reproducibility of genetic variations detection between MIPs-NGS and TSCA-NGS. We considered only the regions defined in Table in S1 Table as target regions. Regions with high variation rates in exons or mapping at the ends of the reads were excluded to avoid mismatches due to the low quality of reads (Table 1). Variants covered less than 30x, with reads supporting the variant call less than 20%, or that were present in more than 90% of the tested samples were filtered-out to exclude potential artifacts. For the cross-comparative analysis, the overlap between each approach was determined. Disconcordant calls were qualitatively and quantitavely assessed.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Performance comparison between MIPs-NGS and TSCA-NGS.

https://doi.org/10.1371/journal.pone.0238467.t001

Results

MIPs-NGS versus TSCA-NGS workflow and costs

MIPs-NGS and TSCA-NGS showed a high multiplexing level and require low sample DNA input (between 50–250 ng). Both protocols are straightforward, but hands-on time and the methodological complexity should be taken into consideration. MIPs-NGS is easier and less time-consuming than TSCA-NGS. The workflow for TSCA-NGS library preparation pass through 86 steps, combining enzymatic and PCR reactions, purification and normalization steps with filer plates and magnetic beads, while MIPs-NGS only required 7 steps to obtain the final library pool [16, 24, 25]. The hands-on time to create a library was 45 minutes for MIPs-NGS versus 130 minutes for TSCA-NGS, regardless of sample size, DNA purity, and concentration of starting material (Table 2 and data in S2 Fig).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Comparison of recommended DNA input, probe features, sequencing kits, sample processing time for MIPs-NGS and TSCA-NGS.

https://doi.org/10.1371/journal.pone.0238467.t002

In this study, costs per sample, including probes, chemical reagents and efforts was lower for MIPs-NGS compared to TSCA-NGS. However, when a large number of patients, e.g. 500 samples were included, the costs per sample for MIPs-NGS could drop to > 5 times as cheap as TSCA-NGS (Table 3).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 3. Price comparison MIPs-NGS and TSCA-NGS on different sequencing platforms.

https://doi.org/10.1371/journal.pone.0238467.t003

Discussion

Selection of an approach for screening a panel of genes for most laboratories depends on a wide range of criteria, including clinical use of the test, panel size, sensitivity and specificity for the genetic regions of interest, expected number of patients to be tested, turnaround time, approach flexibility and scalability, available equipment, work flow, costs, technical expertise, and availability of bioinformatic support.

In this study, we developed a MIPs-NGS and TSCA-NGS targeted re-sequencing panels for nine sodium channel genes (SCN3A, SCN8A-SCN11A, and SCN1B-SCN4B) known to be associated with neuropathic pain [19, 27], and compared their sensitivity, specificity, targeting efficiency, performance and cost effectiveness. The average targeted regions coverage was 97.3% in MIPs-NGS, and 93.9% in TSCA-NGS. Sensitivity, specificity and performance of both methods were comparable (Table 1). However, MIPs-NGS has a more versatile assay design, and is flexible and cheaper than TSCA-NGS (Tables 2 and 3).

Molecular Inversion Probes (MIPs), also known as padlock probes, belong to the category of molecular techniques that capture sequences by circularization. This technology, first described in 1994, was initially developed for multiplex target discovery and SNP genotyping [28, 29], and has recently been combined with next generation sequencing. TSCA-NGS is an amplicon-based approach for targeted re-sequencing. This approach is based on the design of synthetic oligonucleotides (probes), with complementary sequence to the flanking regions of the target DNA to be sequenced. Amplicon-based sequencing approaches are cheaper than Sanger sequencing, WES and WGS, and characterized by high specificity and deep coverage. Moreover, MIPs-NGS and TSCA-NGS have been successfully employed with good-quality DNA sources such as blood or frozen tissues, and with more challenging samples extracted from formalin-fixed and paraffin-embedded tissues [30].

Development of a sensitive diagnostic custom targeted NGS enrichment kit requires proper design of specific primers or probes for candidate genes. For MIPs-NGS and TSCA-NGS probe design, we used the freely available MIPgen software (http://shendurelab.github.io/MIPGEN/) and Design Studio software of Illumina (https://designstudio.illumina.com/), respectively. Both software tools were user-friendly, and simplify the probe designing process based on the criteria defined by the user.

MIPs-NGS and TSCA-NGS showed high multiplexing level, low DNA input requirements (50-250ng) and no need for DNA shearing compared to other targeted enrichment methods [15, 31].

MIPs-NGS and TSCA-NGS have, in general, straightforward laboratory workflows, however TSCA-NGS pass through many steps compared to MIPs-NGS (data in S2 Fig). The hands-on time per library were 2.9 times shorter for MIPs-NGS compared to TSCA-NGS (Table 2 and data in S2 Fig). In addition, development of automated reaction setups is feasible for MIPs-NGS approach since small number of enzymatic reactions and processing steps are required to achieve targeted region capture and sample barcoding. Such an automated laboratory workflow was recently established for smMIP enrichment to detect genetic variations in BRCA1 and BRCA2 [32]. Furthermore, this study concluded that smMIPs-NGS has a superior accuracy and turnaround time compared to other genetic testing methods for gene panels or targeted regions [32].

The performance and capture efficiency of MIPs-NGS and TSCA-NGS targeted approaches for SCN3A, SCN8A-SCN11A, and SCN1B-SCN4B was high (98.9% for MIPs-NGS versus 96.7% for TSCA-NGS) and variation per read low; both required for a reliable variant detection. No sequence reads were obtained by MIPs-NGS and TSCA-NGS for the first coding exon of SCN1B, SCN4B and SCN11A. These regions have a high GC content (>50%) which influence the capture efficiency. To avoid false positive variant calling, data of seven MIPs (1678 bp) versus three TSCA probes (794 bp) should be excluded from the analysis. To provide a full coverage for SCN3A, SCN8A-SCN11A, and SCN1B-SCN4B, 2321 bp (6.1%) for MIPs-NGS versus 2467 bp (6.4%) for TSCA-NGS should subsequently be analysed by Sanger sequencing (costs per sample excluding VAT are €29.7 for MIPs-NGS versus €24.75 for TSCA-NGS).

The utility of custom targeted enrichment NGS panels may be limited in some contexts because the kits are inflexible in terms of adding or excluding targeted regions and usually expensive [33]. In many diagnostic laboratories, adjusting existing diagnostic panels have gained popularity specially in investigating genetically heterogeneous disease [34]. MIPs-NGS approach offered a higher degree in flexibility and in optimizing the probe performance compared to TSCA-NGS and other targeted enrichment NGS methods [2, 35]. All MIPs were ordered individually, so each individual MIP can be combined into various panels and can be added to an existing pool. This means that changing the MIP kit content as the regions of interest change over time and to keep up with the ever-increasing numbers of diagnostics requests is possible. Moreover, the performance of poorly performing probes can be improved by individual MIP rebalancing or by redesigning new probes. In our study the performance of MIPs-NGS was improved by adding a 5 fold concentration of each poorly performing probe (n = 28 MIPs) to the existing probe pool without the need to order new probes. TSCA-NGS is less flexible because changing the kit content or optimizing specific probe(s) performance is dependent on manipulations by the supplying company. However, in this study probe dosage adjustment was not needed.

Applying WGS as a comprehensive clinical testing approach is at this time still too expensive. Sequencing the entire human genome is associated with the generation of a massive amount of data and the need of complex downstream data analyses and variants interpretation. The majority of known disease-causing variants are located in exons; thus WES and targeted gene panels are often used to identify clinically relevant variants in known disease genes [36]. WES -nearly- includes all the exonic regions of the human genes, so using WES eliminates the need to select the gene panel content, and exonic variants are more comprehensively assessed. However, the amount of data generated by WES is significantly higher than targeted gene panels. In contrast to WES [36], targeted NGS enrichment gene panels are often used for re-sequencing a selective group of phenotype-specific disease genes in common genetic disorders. This results in more manageable and easier data interpretation, usually with higher quality and depth at lower costs [2].

MIPs-NGS has a higher performance and capture efficiency for SCN3A, SCN8A-SCN11A, and SCN1B-SCN4B compared to the commonly used Agilent SureSelectXT Human All Exon v5 (Illumina) WES enrichment kit (Table in S2 Table). Similar to our data (Tables 1 and 3), several studies have demonstrated that MIPs-NGS is an efficient and inexpensive method for detection of genetic variations in multiple diseases [25, 32, 37, 38]. For example, Zhang et al. and Pérez Millán et al. presented in their studies the use of MIPs based NGS panels to detect pathogenic mutations in early-onset colorectal cancer patients and in patients with hypopituitarism, respectively [37, 38]. Another potential application of MIPs-NGS is to target noncoding and low covered regions in WES.

High throughput NGS platforms such as NextSeq 500/550 and NovaSeq6000 (Illumina, Inc., San Diego, CA, USA) have enabled reductions in sequencing costs (see Table 3 for comparison costs of sequencing SCN3A, SCN8A-SCN11A, and SCN1B-SCN4B on MiSeq versus NextSeq 500/550 and NovaSeq) and time that have resulted in wide use in diagnostic laboratories. Nevertheless, full capacity high-throughput sequencing runs must be achieved to get the best cost-efficiency (based on maximum number of available sample indexes, 384 samples/run for MIPs-NGS versus 96 samples/run for TSCA-NGS). MIPs-NGS and TSCA-NGS can be easily customized to individual needs. Multiple different targeted MIPs-NGS can be combined cost efficient in one Mid—(48–96 samples) or High output run (96–384 samples) without affecting data quality and coverage for each individual sample by which in diagnostics a turnaround time of 2–4 weeks/sample from arrival to reporting can be achieved.

When designing specific targeted enrichment NGS experiments, total number of samples involved, and cost of probes and reagents must be considered. Our data show that the costs of MIPs-NGS are lower than TSCA-NGS. Using MIPs with a gap-fill length of 220–230 nt [24] and TSCA probes with a region of interest of 425 nt, the cost of MIPs-NGS compared to TSCA-NGS is around 5 times lower per sample (€20.7 for MIPs-NGS versus €113.4 for TSCA-NGS) when 500 samples and more are planned to be tested with MIPs-NGS (Table 3).

In conclusion, our results provide a validation and performance assessment of MIPs-NGS and TSCA-NGS as reliable methods for variant detection in a disease-specific subset of genes. Our results suggests that MIP-NGS is a more flexible and less expensive method for detection of genetic variations and is a reliable screening approach for laboratories involved in diagnostic service for patients with pain-related disorders.

Supporting information

Acknowledgments

Lead author PROPANE Study Group: G. Lauria, Neuroalgology Units, Fondazione IRCCS Istituto Neurologico “Carlo Besta” Milan, Italy and Department of Biomedical and Clinical Sciences "Luigi Sacco", University of Milan, Milan, Italy, Email: Giuseppe.LauriaPinter@istituto-besta.it.

Members Propane Study Group: Bianca M. de Greef (Maastricht UMC+, Maastricht, the Netherlands); Jo M. Vanoevelen (Maastricht UMC+, Maastricht, the Netherlands); Ivo Eijkenboom (Maastricht University, Maastricht, the Netherlands); Michela M. Taiana (IRCCS Foundation Neurological Institute “Carlo Besta”, Milan, Italy); Raffaella Lombardi ((IRCCS Foundation Neurological Institute “Carlo Besta”, Milan, Italy); Daniele Cazzato (IRCCS Foundation Neurological Institute “Carlo Besta”, Milan, Italy); Andrea Zauli (San Raffaele Scientific Institute, Milan, Italy); Ferdinando Clarelli (San Raffaele Scientific Institute, Milan, Italy); Ignazio D. Lopez (San Raffaele Scientific Institute, Milan, Italy); Angelo Quattrini (San Raffaele Scientific Institute, Milan, Italy); Mitra Tavakoli (University of Manchester, Manchester, United Kingdom); Dimos Kapetis (IRCCS Foundation Neurological Institute “Carlo Besta”, Milan, Italy); Markos N. Xenakis (Maastricht University, Maastricht, the Netherlands); Massimo Mantegazza (University of Côte d’Azur, Sophia Antipolis, France); Francesco Battiato (German Diabetes Center, Düsseldorf, Germany); Alexander Strom (German Diabetes Center, Düsseldorf, Germany); Sandrine Cestèle (University of Côte d’Azur, Sophia Antipolis, France); Oana Chever (Collège de France, Paris, France).