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

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.

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).

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.

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.

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.