MotorPlex provides accurate variant detection across large muscle genes both in single myopathic patients and in pools of DNA samples

Abstract

Mutations in ~100 genes cause muscle diseases with complex and often unexplained genotype/phenotype correlations. Next-generation sequencing studies identify a greater-than-expected number of genetic variations in the human genome. This suggests that existing clinical monogenic testing systematically miss very relevant information.

We have created a core panel of genes that cause all known forms of nonsyndromic muscle disorders (MotorPlex). It comprises 93 loci, among which are the largest and most complex human genes, such as TTN, RYR1, NEB and DMD. MotorPlex captures at least 99.2% of 2,544 exons with a very accurate and uniform coverage. This quality is highlighted by the discovery of 20-30% more variations in comparison with whole exome sequencing. The coverage homogeneity has also made feasible to apply a cost-effective pooled sequencing strategy while maintaining optimal sensitivity and specificity.

We studied 177 unresolved cases of myopathies for which the best candidate genes were previously excluded. We have identified known pathogenic variants in 52 patients and potential causative ones in further 56 patients. We have also discovered 23 patients showing multiple true disease-associated variants suggesting complex inheritance. Moreover, we frequently detected other nonsynonymous variants of unknown significance in the largest muscle genes. Cost-effective combinatorial pools of DNA samples were similarly accurate (97-99%).

MotorPlex is a very robust platform that overcomes for power, costs, speed, sensitivity and specificity the gene-by-gene strategy. The applicability of pooling makes this tool affordable for the screening of genetic variability of muscle genes also in a larger population. We consider that our strategy can have much broader applications.

Introduction

Muscle genetic disorders comprise about 100 different genetic conditions [1],[2], characterized by a clinical, genetic and biochemical heterogeneity. The molecular diagnosis for myopathic patients is crucial for genetic counseling, for prognosis and for available and forthcoming mutation-specific treatments [3]-[5]. In addition, patients that share the same mutation may have a different type of muscle affection with the selective involvement of other muscle compartments or myocardial damage. Thus, the primary defect may be modified or not by additional and variable elements that may be genetic or not. The most severe cases of congenital or childhood-onset myopathies often result from mutations in genes encoding proteins belonging to common pathways [6]. To provide a clue to address genetic testing, a muscle biopsy is often required that may be useful, but not well accepted by patients. The single gene testing can be diagnostic only in patients with most recognizable disorders. In unspecific cases of muscular diseases, however, no effective methodology has been developed for the parallel testing of all disease genes identified so far [7].

Next-generation sequencing (NGS) is changing our view of biology and medicine allowing the large-scale calling of small variations in DNA sequences [8]. In the last few years, the whole-exome sequencing (WES) and whole-genome sequencing (WGS) have received widespread recognition as universal tests for the discovery of novel causes of Mendelian disorders in families [9]. The power to discover a novel Mendelian condition increases with the family size, even if successful studies, identifying novel disease genes from multiple small families with the same phenotype, have been published [10]. Structural and copy number variations are not well detected by NGS technologies [11]-[14]. However, the WES/WGS use for the clinical testing of isolated cases is still debated. First, there are ethical issues linked to the management of the incidental findings [15]. The second limitation is given by the practical problem that the coverage is usually too low for clinical diagnosis. Hence the cost-effectiveness is reduced, considering that WES/WGS may require either numerous validation procedures, mainly based on conventional PCR and Sanger sequencing reactions [16]. Innovative strategies of clinical exome sequencing at high coverage have been described [17], but the cost for a single patient is still too high for routine diagnosis. Thus, there is still space for targeted strategies [18] and the HaloPlex Target Enrichment System [19] represents an innovative technology for targeting, since it uses a combination of eight different enzyme restriction followed by probe capture. It permits a single-tube target amplification and one can accurately predict the precise sequence coverage in advance. We have developed a NGS targeting workflow as a single testing methodology for the diagnosis of genetic myopathies that we named Motorplex. Here we demonstrate the high sensitivity and specificity of Motorplex. We challenged our platform against complex DNA pools. Even with this complexity, Motorplex kept producing reliable data with high sensitivity and specificity values. Furthermore, pooling reduced the cost of the entire analysis at negligible values, implementing applications for large studies of populations [16],[20].

Materials and methods

Patients

Encrypted DNA samples from patients with clinical diagnosis of nonspecific myopathies, congenital myopathy, proximal muscle weakness or limb-girdle muscular dystrophy (LGMD) were included. The Italian Networks of Congenital Myopathies (coordinated by C.B. and F.M.S.) of LGMD (by F.M. and G.P.C.) were involved together with a large number of other single clinical centers. We asked all them the possibility to share more clinical and laboratory findings, when necessary. We also requested to provide information on familial segregation and previous negative genetic tests. Internal patients signed a written informed consent, according to the guidelines of Telethon Italy and approved by the Ethics Committee of the “Seconda Università degli Studi di Napoli”, Naples, Italy.

DNA samples were extracted using standard procedures. DNA quality and quantity were assessed using both spectrophotometric (Nanodrop ND 1000, Thermo Scientific Inc., Rockford, IL, USA) and fluorometry-based (Qubit 2.0 Fluorometer, Life Technologies, Carlsbad, CA, USA) methods.

In silico design of MotorPlex

We included in the design all the 93 genes that are universally considered as genetic causes of nonsyndromic myopathies (Additional file 1: Table S1). In particular, we only selected genes determining a primary skeletal muscle disease, such as underlying muscular dystrophies, congenital myopathies, metabolic myopathies, congenital muscular dystrophies, Emery-Dreifuss muscular dystrophy, etc. We therefore excluded loci associated with other neuromuscular and neurological disorders such as congenital myasthenias, myotonic dystrophy, spinal muscular atrophy, ataxias, neuropathies, or paraplegias for which differential diagnosis may be clinically possible. For each locus, all predicted exons and at least ten flanking nucleotides were always included in the electronic design by the custom NGS Agilent SureDesign webtool. Setting the sequence length at 100×2 nucleotides, the predicted target size amounted to 2,544 regions and 493.598kb. Around 20% of the target is represented by TTN coding regions.

NGS workflow

For library preparation of single samples, we followed the manufacturer’s instructions (HaloPlex Target Enrichment System For Illumina Sequencing, Protocol version D, August 2012, Agilent Technologies, Santa Clara, CA, USA). We started using 200ng of genomic DNA and strictly followed the protocol, with the exception that restricted fragments were hybridized for at least 16–24 hours to the specific probes. After the capture of biotinylated target DNA, using streptavidin beads, nicks in the circularized fragments were closed by a ligase. Finally, the captured target DNA was eluted by NaOH and amplified by PCR. Amplified target molecules were purified using Agencourt AMPure XP beads (Beckman Coulter Genomics, Bernried am Starnberger See, Germany).

The enriched target DNA in each library sample was validated and quantified by microfluidics analysis using the Bioanalyzer High Sensitivity DNA Assay kit (Agilent Technologies) and the 2100 Bioanalyzer with the 2100 Expert Software. Usually 20 individual samples were run in a single lane (250M reads), generating 100-bp paired end reads.

For Pool-Seq experiments, equimolar pools of 5 or 16 DNA samples (detector and scouting pools) were created and 200ng of each pool was used for the HaloPlex enrichment strategy. Sixteen detector and five scouting pools were usually run in a single HiSeq1000 lane.

Targeted sequencing analysis

The libraries were sequenced using the HiSeq1000 system (Illumina inc., San Diego, CA, USA). The generated sequences were analyzed using an in-house pipeline designed to automate the analysis workflow, composed by modules performing every step using the appropriate tools available to the scientific community or developed in-house [21]. Paired sequencing reads were aligned to the reference genome (UCSC, hg19 build) using BWA [22], sorted with Picard (http://picard.sourceforge.net) and locally realigned around insertions-deletions with Genome Analysis Toolkit (GATK) [23]. The UnifiedGenotyper algorithm of GATK was used for SNV and small insertions-deletions (ins-del) calling, with parameters adapted to the Haloplex-generated sequences. The analysis of pools was performed with UnifiedGenotyper as well, adapting the ploidy parameter to the number of chromosomes present in the samples (10 for the detector and 32 for the scout pools) and the minimal ins-del fraction parameter accordingly. The called SNV and ins-del variants produced with both platforms were annotated using ANNOVAR [24] with: the relative position in genes using RefSeq [25] gene model, amino acid change, presence in dbSNP v137 [26], frequency in NHLBI Exome Variant Server (http://evs.gs.washington.edu/EVS) and 1000 genomes large scale projects [27], multiple cross-species conservation [28],[29] and prediction scores of damaging on protein activity [30]-[33]. The annotated variants were then imported into the internal variation database, which stores all the variations found in the re-sequencing projects performed so far in our institute. The database was then queried to generate the filtered list of variations and the internal database frequency in samples with unrelated phenotype was used as further annotation and filtering criteria. The alignments at candidate positions were visually inspected using the Integrative genomics viewer (IGV) [34]. We selected from the database the non-synonymous SNVs and ins-del, with a frequency lower than 2%, which was followed by manual inspection and further filtering criteria based on the presence in unrelated samples of the database, on the presence in the other samples of the Motorplex experiment and on the conservation of the mutations, with a final selection of rare, possibly causative, variations per individual.

Results

Validation study of MotorPlex

To design MotorPlex we used a straightforward procedure. Briefly, disease genes causing a muscular phenotype, including the biggest genes of the human genome, like titin (TTN) or dystrophin (DMD), were selected. The target sequences, corresponding to 0.5Mbp were enriched by the HaloPlex system (see Materials and methods). To validate MotorPlex, we created a training set of twenty DNA samples belonging to patients (15 males and 5 females) affected by different forms of limb-girdle muscular dystrophy or congenital myopathy (Additional file 2: Table S2) and compared with data from whole exome sequencing (WES) (Figure 1). For each sample, about 98% of reads generated (Figure 1a and Additional file 3: Table S4) were on target (compared to 88% obtained by WES) and fewer than 0.5% of targeted regions were not covered (about 15% of human exons are not analyzed by WES, Additional file 4: Figure S1). Moreover, more than 95% of targeted nucleotides were read at a 100× depth and a 500× depth was obtained for 80% of these; on the contrary, by performing a WES analysis, fewer than 70% of exons were covered at 20× (Figure 1b). From previous amplicon Sanger sequencing from these samples, we knew about 84 variants in 17 different genes (Additional file 5: Table S3). All these known variants were correctly called and no additional change was seen within the sequenced target (100% sensitivity and specificity). Moreover, to assess the reproducibility of the targeted enrichment and the subsequent NGS workflow, the same sample (43U) was analyzed twice. After filtering, variants were always confirmed, including the putative causative one (Table 1). Outside the Sanger coverage, 4,991 additional variations were called (Additional file 6: Table S5).

Figure 1

A comparison between MotorPlex and a Whole Exome strategy (WES) demonstrates the better performance of the targeted strategy. (a) 97.75% of reads generated in a MotorPlex experiment fall in the regions of interest and only 0.67% of targeted regions are not sequenced. On the contrary, for WES 88.66% of reads are on target and 14.89% of targeted exons are not effectively covered. (b) The percentage of targeted regions covered at high depth by MotorPlex is higher than that obtained by WES. In particular, 96.01% and 81.6% of regions are, respectively, covered at 100x and 200x by using MotorPlex versus 35.49% and 1.90% by WES.

Table 1 List of pathogenic variants

Validation study of double-check pooling

To challenge MotorPlex to be applied to large studies on thousands of patients and/or to detect mosaic mutations, we designed a combinatorial pooling strategy. After some initial attempts with pools of identical sizes, we changed our strategy. The general arrangement was to have the same sample in two different independent pools, composed of two exclusive combinations of samples (Figure 2). This permitted us to identify both the rare variations and the sample mutated. In particular, the pools were organized in two groups: the “detector pool” only containing five samples (10 alleles) that had the purpose of detecting variations with the optimal sensitivity and the “scout pool” composed of 16 samples (32 alleles) that confirmed the variation(s) and attribute them univocally to distinct DNA samples (Additional file 7: Figure S2; Additional file 8, Table S6). We paid attention each time to include the index cases alone, excluding related family members.

Figure 2

NGS targeting workflow. Ninety-three disease genes causing a muscular phenotype were selected. To cover all their exons and the ten flanking bases, an enrichment strategy, based on HaloPlex system, was designed. DNA samples of 80 patients were analyzed twice in an independent manner, using a combinatorial pooling scheme. As requested by HaloPlex protocol, DNA samples were digested, barcoded and amplified. The 80 samples were run at the same time in a single lane of the flow cell of HiSeq 1000. The following data analysis allowed us to detect putative causative variants validated by Sanger sequencing.

To validate this arrangement, we selected five samples that we previously sequenced individually and called 1,235 variations. We pooled them in the same detector pool (P9) and then reanalyzed in different scout pools. Impressively, in pool P9 we called 1,232/1,235 variations belonging to the individual samples, calculating the sensitivity value at 99.8%. The three missing variations (an insertion in RRM2B and two point variants in TTN) were located in regions with lower coverage. On the contrary, no variation was called in pool 9 in addition to those of individual samples, demonstrating the absence of false positives and artefacts due to the pooling strategy. Another two samples from the training set were inserted in another two detector pools, showing similar results.

We then confirmed 223/230 (97%) variations tested by Sanger sequencing, thus providing the specificity value of the method. Moreover, the combined use of detector pools and scout pools allowed us to “clean” the results. 50% of off target variations (n=1,291), in fact, were not called in the scout pools and were easily filtered off during bionformatic analysis. In addition, about 25% of variants in low covered regions (<500 total reads), representing in a large percentage false positive calls, were similarly filtered off because they were not detected in the scout pools (Additional file 9: Figure S3).

Variants and interpretation

The targeted analysis of 93 genes showed a total of 23,109 rare variants (<0.01 frequency) in 173 patients (1.4 variants/gene/patient). To provide a preliminary interpretation in relationship with the clinical suspicion, we set bioinformatic filters that weigh the variant class (missense, indel, stopgain or stoploss), the calculated frequency in public and internal databases and the annotation as causative variants. Finally, we reconsidered critically the correspondence with the clinical presentation, the age at onset and the segregation in familial cases.

In detail, we identified 52 patients (52/177=29%) with variants of likely pathogenicity or predicted to affect function (Table 1 and Additional file 10): most of them (38/52=73%) had known or truncating variants (indel, stopgtain or stoploss). Five patients (5/52=9.6%) showed a novel variant in addition to a pathogenic allele in a recessive gene. The remaining samples (9/52=17%) had novel variants that are predicted to affect function in genes fitting with the clinical suspicion.

In other 56 samples (56/177=32%), we identified potential causative variants (Table 2 and Additional file 10). In these cases, there was only a partial correspondence with the clinical phenotype. For example, a number of variants had been previously associated with cardiomyopathy, but their pathogenic role in congenital myopathy or in LGMDs was not yet established. To the group belong patients having two rare variants in TTN gene or at least one variant in COL6A1, COL6A2, COL6A3, SYNE1, SYNE2 and FLNC genes. These molecular findings in these 56 samples were not considered strictly disease-causing and further tests are required.

Table 2 Variants of unknown significance (Vous)

The most surprising finding was, however, the presence of additional damaging or potential damaging variants in 16 patients of the first two groups (23/108=21%) in whom other pathogenic variants or variants of uncertain significance had already been identified. These variants, if they had been detected alone in the context of a single gene testing, would have been considered as causative.

The third group includes 26 patients (26/177=15%) in which we discovered a single truncating variant (or a known disease-associated variant) in a recessive gene that is compatible with the phenotype. The second allele may carry a RNA splicing defect that is generally not predictable by DNA sequencing or, also, a variation in not investigated promoters or regulatory regions.

Discussion

In the last decade, a remarkable progress has been made in discovering new disease genes and differentiating similar muscle disorders [1],[2]. This growing genetic heterogeneity highlights the problem of a very complex diagnosis [35]. Furthermore, genome sequencing studies suggest that the clinical genetic test may be incomplete not only when the causative mutation is missing, but also when the genotype/phenotype correlation appears weak. This is particularly true when the familial recurrence is unclear, with some relatives that only share minor affections. In families with patients who are more severely affected, this “grey area” is problematic for both genetic counselling and forthcoming mutation-specific treatments. However, this represents the proper challenge for the new genomic, high-throughput technologies: the power of discovery has been dramatically boosted by the introduction of the next-generation sequencing (NGS) techniques [13],[36]-[38]. In the NGS era, the genetic testing is going to move from few candidate genes to broader panels of genes [39] or, ultimately, to the entire genome. This will have consequences on the diagnostic flowchart: NGS tests may represent the first tier test, preceding biopsy and other invasive procedures.

We have applied both WES and targeted approaches to the diagnosis of genetic disorders of muscle and collected DNA samples of patients without diagnosis and realized that NGS technology can be helpful for clinical diagnostics, provided that a suitable tool is created. We traced an ideal profile of it. This tool should fulfil the following requirements [16],[20]: 1) to be cost-effective and thus applicable to a large number of patients and normal individuals, 2) to be robust in the terms of target reproducibility, 3) to be specific and sensitive with a limited need for further validation steps, 4) to be large enough to include all relevant genes and, finally, 5) to be easily upgradable in view of novel discoveries. Here we demonstrate the ability to generate this complex targeting and to fulfil all these requirements. We decided to use Haloplex as the enrichment technology. Haloplex first digests DNA using eight different combinations of endonucleases. Our experience suggests that this approach is more reproducible and accurate than the random mechanical DNA fragmentation. In addition, the capture is independent of the target base composition and is predictable from the probe design phase. As a proof of specificity and efficiency, we show that less than 2% of reads generated by Motorplex are off-target, in comparison with >12% of WES. This factor further improves the cost-effectiveness of the approach. This platform, based on eight different digestions and hybridization, is more accurate, reproducible and sensitive in comparison with other published methods [34]. We have designed the MotorPlex to detect variations in 93 muscle-disease genes and assayed 177 pre-screened DNA samples from myopathic patients. It is important to consider that these are all patients with zero mutations so far detected, even if most of them have been lengthily studied using a gene-by-gene sequencing approach. The high coverage and depth obtained permitted us to detect variations in most genes with sensitivity comparable with Sanger sequencing. According to our conservative NGS data interpretation, in 52 patients (29%) the diagnosis is complete. However, the detection rate will grow after a further molecular characterization of putative pathogenic variations in a second group of 56 patients. In addition, there are 26 samples (15%) that have defects in one single allele associated with a recessive condition. We predict that most of these can carry an elusive hit on the other allele such as splicing defects or copy number mutation(s). A percentage of 15%, in fact, is a usual value for disease-causing variants not detectable by sequencing.

The most interesting and quite surprising finding is, however, the very high number of rare damaging variants identified and first the cases (26/177) with more damaging variants in other genes in addition to those classified as causative. These additional variants may have a potential modifier effect. This percentage of these genetically complex patients may be higher, if we consider that many other important muscular genes (even if not disease-causing) can also carry damaging alleles. We can easily predict that a broader NGS approach could strengthen this observation. We hypothesize that the intrafamilial and interfamilial phenotypic differences may be frequently related to the combinations of multiple disease-causing alleles, more than to SNPs or CNVs. The so-called “modifier gene variants” could be individually rare, but collectively common. A comprehensive view of all the genes involved in a pathological process helps to point out these alleles having a minor but probably not negligible role in the disease aetiology.

The ultimate goal of MotorPlex is given by the pooling performances. The specificity and sensitivity values are very high and quite similar to those obtained in singleton testing and, above all, the diagnostic rate is not affected. The potential applications of pooling are just in large studies of complex and non-Mendelian disorders when a large number of samples have to be analyzed to improve the statistical power [40]. Considering our finding of multiple damaging variants in disease genes, these large studies are just around the corner. In addition, MotorPlex may discover low-allelic fraction variants in single samples, as in somatic mosaicisms. The pooled MotorPlex is likewise the cheapest genetic test (Table 3) ever presented that is able to screen 93 complex conditions at the cost of a few PCR reactions.

Table 3 Predicted enrichment costs and workload for single and pooled DNA samples

In conclusion, we here demonstrate that MotorPlex can be used to identify accurately all DNA variants also in huge muscle genes: the platform overcomes for sensitivity and coverage the WES approach. In addition, Pool-Seq may be the first option to perform cost-effective population studies to understand polygenic conditions. We think that similar protocols could be designed to extend the NGS applications to other studies for human genetics, as well as for disease prevention, nutrition, forensics and many others.