Sufficient coverage of PKD1 and other genes known for cystic and polycystic kidneys as a prerequisite for diagnostic testing

In total, we sequenced 55 samples, 10 samples on the MiSeq and 45 samples on two lanes of a HiSeq 1500 system with 2 × 150 bp paired-end runs. For the whole target region including genes for cystic and polycystic kidney disease and potential phenocopies (in total, n = 40) (S1 Table) enriched and sequenced in parallel to PKD1 a mean target coverage of 650x was achieved on the MiSeq and of 1443x on the HiSeq system on average across all samples with about 96.6% and 97.4% of the target regions covered with at least 20x, respectively (S2 Table). The average coverage of the completely targeted genomic region of the PKD1 gene was 914x (91.8% >20x) on the MiSeq and 1993x (92.9% >20x) on the HiSeq (Fig. 1A) with gaps in the intronic sections. The 46 coding exons of PKD1 were covered with an average coverage of 1198x and 2584x, respectively, nearly covering the complete coding sequence with 20x (98.6 and 99.1%). Moreover, nearly all coding positions in all exons were covered far more than 100x in all samples (Fig. 1B). As expected due to major hassles also known from conventional Sanger sequencing, the only region that could not be completely sequenced was the first half of exon 1 due to its high GC content and despite thorough rebalancing with 50.6% and 23.5% on average being 20x covered on the HiSeq and MiSeq, respectively. Nearly all target genes beside PKD1 were on average covered >20x for 100% of the coding region with only minor exceptions due technical limitations (S3 Table).

SNV and small indel detection by our capture-based NGS approach

NGS data of 53 PKD1 SNV (single nucleotide variant)/ indel (insertions and deletions) positive patient samples were compared to results of previous analyses of these patients by conventional LR-PCR followed by Sanger sequencing. All but one of the 160 different variants (681 of 683 total variants, Table 1) virtually distributed evenly over all 46 PKD1 exons (Fig. 1A) and including different mutation types (various SNVs, short insertions up to a 15 bp-duplication as well as short deletions with a maximal 12 bp-deletion) were identified (S4 Table) mainly with standard settings. All definitely pathogenic mutations, in total 31 nonsense mutations and indels, were clearly detected. One variant not detected by MiSeq sequencing located in exon 1 that was found to be insufficiently covered (S5 Table). However, this variant might have been detected when sequencing the sample on the HiSeq instrument (Fig. 1). Among putatively pathogenic variants, which were either predicted pathogenic (S6 Table) or had been reported with indefinite nature or were private to the respective patient, only one could not be detected (S5 Table). This specific indeterminate variant, c.2180T>C (p.Leu727Pro), in exon 11 could in one case not be detected at all due to missing the stringent and invariable filter criteria in variant calling (patient 21) and in a second case only after reducing the detection threshold in a second step analysis. In total, SNVs at only six distinct sites (Fig. 2A) in five exons of the duplicated region were either detected below the threshold of 20% reads at the altered position or were called with different zygosity (S5 Table). Inspection of the alignment showed two different cases: First, in most of these cases variants at these sites exist in at least one pseudogene copy with the complication that reads carrying the exchange map to 100% to a homologous region. Only few read pairs are present in the alignment which map to the master gene as discrimination is enabled by the second read of this pair (S1 Fig.). However, a proportion of all pairs show low mapping quality due to little sequence divergence between master gene and pseudogenes. In the second case two of these critical sites show sequence identity between PKD1 and the pseudogenes at the site of interest. Filtering of reads against their mapping quality resulted in better results only at these two sites maybe due to filtering for pairs carrying the variant, for which mapping discrimination resulted in higher mapping quality (S2 Fig.). However, in general, the challenge of variant detection at these few sites could be coped with a second-step analysis by simple adjustment of the detection threshold to 8% alternative reads in all but the two above discussed false negative sites. It should be emphasized that relaxed filter criteria did not result in an increased detection of false-positive calls showing high specificity of the applied mapping algorithm even at a lower stringency. Interestingly, the 17 allegedly false positive calls not been detected by conventional testing restricted to only two sites in exons 10 and 11. Because of sequence identity of the genuine gene with the duplicated regions at these sites any spurious alignment of pseudogene regions causing their calls may be ruled out. As these two false-positive variants were detected in NGS only slightly above the threshold of 20% reads at the altered position and as they were not detectable by Sanger sequencing even by using alternative primer pairs, they most probably represent recurring sequencing artefacts specific to the whole workflow, which are to be filtered out by their internally measured allele frequency. Taken together, the approach showed a very high sensitivity (99.7%) and specificity (99.8%) (S4 Table) comparable to conventional testing methods.

Variant simulation in the duplicated PKD1 region with a high variant density underscores sensitivity of our approach

In our cohort almost all exons of the PKD1 gene were covered by real variants. To have a closer look at gaps that escaped evaluation and to rule out further coding segments in the duplicated region, in which detection of variants might be challenging due to the high homology between PKD1 and its pseudogenes at and around a site of interest, SNVs and indels in the PKD1 gene were simulated using Wgsim by having at least one alteration within each 100 bp of coding sequence in the duplicated region (exons 1–33). Simulated reads were generated with a coverage of 1000x and combined with unmodified reference genomic sequence of all six pseudogenes (PKD1P1–6) with equal coverage. The combined reads were finally mapped and analysed using our in-depth bioinformatic pipeline (S3 Fig.). Thereby, a total of 241 different mutations were simulated with an even distribution over all 33 exons of the duplicated region (Fig. 2B, S7 Table), which is an ever higher variant density than for the 126 different changes from the validation cohort in this region. All except one mutation in exon 5 could be detected at standard filter criteria (20% alternative reads) resulting in 99.6% sensitivity without adjustment of the detection threshold. Interestingly, the simulated nucleotide exchange at position c.1155 detected only in 15.5% of all reads is located close to the most common critical real variant at position c.1119 discussed above demonstrating that variant detection at this segment is challenging due to competitive mapping with pseudogene regions emphasizing how realistic our simulation reflected the real world scenario.

Potential CNV detection in all target genes

We utilized VarScan to evaluate the potential of our approach for CNV analysis. Two patients harbouring previously known deletions in PKD1 (patient 30 carrying a single-exon deletion of exon 22 and patient 46 bearing a genomic deletion comprising exons 15–21) were used as positive controls. Six samples that were previously proven by MLPA to have normal copy numbers for PKD1 served as negative controls. As expected from our results for SNV calling, regions with evaluated competitive mapping in the duplicated region (five exons) were challenging to analyse and generated false-positive CNV calls in these exons. Therefore these and recurrent CNV artefacts as for exons 1 and 42 attributable to insufficient or varying coverage of these regions well-known for their complexity were excluded from the analysis. Applying those stringent filter criteria, a deletion of exon 15 in patient 46 and a deletion of exon 22 in patient 30 could be clearly detected, showing the CNV convincingly against all control samples (Fig. 3 and S4 Fig., S8 Table). When lowering the detection threshold to against 50% control samples and with the anchor point of exon 15 being most probably deleted, the extent of the above multi-exon deletion encompassing exons 15–21 could be precisely predicted from NGS CNV data. Interestingly, combining this depth of coverage detection method with a split read approach using ClipCrop, 332 reads were obtained which supported a breakpoint in exon 15 (Fig. 3). A Smith-Waterman alignment of the unmapped part of the spilt reads against the entire DNA sequence of PKD1 found the best match in intron 21 known to be repetitive.

The feasibility of the sequence capture approach for CNV analysis in the other targeted genes could be demonstrated by detection of a known heterozygous deletion of the complete NPHP1 gene in patient 3 (Fig. 4D). As seen in previous gene-panel sequencing setups detection of deletions in the HNF1B gene, which might account for up to 50% of HNF1B-mutation positive cases, is also accessible by a sequence capture based NGS approach (S5 Fig.). Overall, CNV detection in the other targeted genes might be even more solid than for PKD1. Further optimization is reasonable before conventional MLPA analysis for PKD1 can be ultimately regarded superfluous in mutation-negative ADPKD patients. However, under stringent filter criteria for PKD1 analysis CNVs present in our patient cohort were correctly detected.

Evaluation of variants in other genes for cystic and polycystic kidney disease

As seen for CNV analysis and as already proven by previous gene-panel based NGS studies parallel analysis of all other targeted genes does not constitute a major technical issue [14],[15],[16]. As proof-of-principle NGS data of our patient cohort were comprehensively analysed for pathogenic aberrations in other genes for cystic and polycystic kidney disease which were all sufficiently covered for diagnostic testing (Table 2). For example, in patient 10 the detection of a mutation in PKD2 (Fig. 4A) underscored the feasibility of the approach for parallel analysis of the second known gene for ADPKD. In this patient with suspected ARPKD and negative family history with normal parental ultrasounds, trans-heterozygosity of a hypomorphic PKD1 allele and a PKD2 mutation might be speculated on, but this hypothesis needs further validation before firm conclusions can be drawn. In patient 5 with unremarkable family history and a phenotype in between ARPKD and ADPKD, a nonsense mutation in HNF1B was identified in trans to a likely hypomorphic PKD1 variant (Fig. 4B). In a number of other patients with early and severe PKD and negative family history suspicious for ARPKD, conventional approaches would have started with Sanger sequencing of PKHD1 which is already laborious and costly due to the huge size of this gene. In some of these patients the identification of a mutation in PKHD1 would have potentially led to the wrong assumption that PKHD1 might be the major disease locus and that genetic testing of other genes is thus dispensable (Fig. 4C). However, only knowledge of the mutational status of other genes such as PKD1 revealed that the main disease locus is a different one in these individuals (S6 Fig.). Conclusively, our data demonstrate the power of the setup for parallel analysis of all genes known to date for cystic and polycystic kidney disease and related disorders in a single step.

Ethics Statement

All samples were obtained with written informed consent. Clinical investigations were conducted according to the Declaration of Helsinki. The study was approved by the institutional review board of the Ethics Committee of the University Freiburg Medical Center.

Samples and preceded Sanger sequencing

DNA samples of 55 patients were included as positive controls in this study, 53 of them were previously sequenced by conventional LR-PCR based Sanger sequencing. Two samples harboured a large deletion in the PKD1 gene previously detected by MLPA (multiplex ligation dependent probe amplification). Genomic DNA was isolated from EDTA blood following standard protocols. Preceded Sanger sequencing had been conducted on long-range PCR (LR-PCR) amplicons as described elsewhere [8,9]. For reassessment of the exons 10 and 11 of PKD1 alternative primer pairs were used as described in [19] (S1 Materials And Methods). Identified sequence variants had been annotated according to the guidelines published by the Human Genome Variation Society.

NimbleGen SeqCap EZ choice library design

The complete genomic region of the PKD1 gene with 1 kb of flanking sequence (NC_000016.9; chr16:2137709–2186899) was targeted by a custom SeqCap EZ choice library (NimbleGen, Madison, Wisconsin, USA). In order to capture the exons in the duplicated region the preferred stringent settings of allowing three close matches in the genome for probe selection were relaxed to up to ten close matches. Due to the high GC content of exon 1 and the low coverage for exon 42 of PKD1 obtained in previous sequence capture approaches we added additional twelve times replicated probes in these regions (rebalancing). In terms of differential diagnosis, the complete genomic sequence of the PKD2 gene and the coding and non-coding exons of 38 other genes for allied human diseases with related clinical manifestations (S1 Table) with 20 bp of flanking genomic sequence were additionally targeted by our design with preferred settings resulting in 820 target regions with 300 kb of target sequence. Rebalancing was also applied to regions showing low depth of coverage in previous panel sequencing (e.g. exon 1 of PKD2). The design covered about 97% of all requested target regions. The coding exons of the PKD1 gene could be completely targeted and about 84% of the whole genomic region was covered. For more details about the NimbleGen design see S1 Materials And Methods. The BED file of the final capture design is available on request.

Sequence capture and next-generation sequencing

Libraries were prepared using NEBNext DNA library prep master mix set for Illumina (New England Biolabs, Ipswich, MA, USA) according to the manufacturer’s protocol. In brief, 1 μg of genomic DNA per sample was sheared to generate DNA fragments of 250–300 bp using the Covaris S2 AFA system (Covaris Inc. Woburn, MA, USA) and ligated to Illumina specific adaptors for multiplexing. Pre-capture amplified samples were pooled—10 samples for MiSeq and up to 23 samples for HiSeq 1500 sequencing—and hybridized to the customized in-solution capture library for 72 hours, subsequently eluted and post-capture amplified by ligation mediated (LM-) PCR. This amplified enriched DNA was used as input for direct cluster generation and sequencing on an Illumina MiSeq system (2 × 150 bp PE reads) or for cluster generation on a cBot system and sequencing in rapid mode (23 samples per lane, altogether up to 46 samples) on an Illumina HiSeq 1500 instrument (2 × 150 bp PE reads) (Illumina, San Diego, CA, USA). A more detailed protocol is available in online S1 Materials And Methods.

Bioinformatic data analysis

NGS data were analysed with an in-house bioinformatic pipeline as previously described [11],[16],[29]. Demultiplexed reads from Illumina MiSeq and HiSeq sequencers were mapped against the hg19 human reference genome using BWA v0.7.8 [30] with the BWA-MEM alignment algorithm and the recommended standard settings and preprocessed with SAMtools v0.1.19 [31] to convert the SAM files into BAM files and sorting them on coordinate order. Duplicate reads were marked with MarkDuplicates from Picard v1.112 (http://picard.sourceforge.net), and the tools RealignerTargetCreator, IndelRealigner, BaseRecalibrator and PrintReads from GATK (Genome Analysis Toolkit) v3.1 software package [32] were applied for a local realignment and base quality score recalibration of the mapped reads. All tools were used with the recommended standard settings. This workflow is in accordance with the best practices from the Broad Institute. Variants were called with the tool UnifiedGenotyper from GATK. Mapping and coverage statistics for the PKD1 regions were calculated using the tool DepthOfCoverage from GATK, as well as the tool CalculateHsMetrices from Picard. The identified variants were annotated using ANNOVAR and the RefSeq gene-based annotation method. After this, all variants were checks against the population frequency databases 1000 Genomes Project (1000g2012apr), Exome Sequencing Project (esp6500si_all). Variants have already been annotated with the dbSNP build 138 database during variant detection with UnifiedGenotyper in a previous workflow step. A more detailed description of the bioinformatics workflow and used scripts are supplied in online S1 Materials And Methods. Annotation (dbNSFP, HGMD), functional prediction and classification of identified variants was conducted as described previously [11],[29]. Additionally, PKD1 variants were annotated with the entry from the ADPKD Mutation Database Version 3.0 (http://pkdb.mayo.edu/). Results from NGS were compared to results from preceded Sanger sequencing (gold standard) in exonic regions ± 30 bp of flanking intronic sequence using bioinformatic scripts and by visual comparison.

For detection of variants in the PKD1 gene the following criteria were applied:

1. Variants were included in the analysis when ≥ 20% of total reads at the position showed the alteration. A variant was called homozygous when ≥ 85% of all reads had the variation.
2. Variants previously identified by Sanger sequencing or annotated in the ADPKD database which could only be detected below the threshold of 20% of total reads were used as indicators for exonic areas in the duplicated region (exon 1–33) where discriminative mapping of reads between master gene and pseudogene was insufficient due to high sequence homology in these reads or where the alternative allele might consequently be underrepresented (S5 Table). Such cases could be evaluated in five exons (5, 11, 17, 21, 26; Fig. 1). If filtering against variants with reads ≥ 20% did not result in detection of a pathogenic alteration, in a second step analysis criteria in the duplicated region of PKD1 (exon 1–33) were relaxed to ≥ 8% to detect candidate pathogenic variants for subsequent confirmation by Sanger sequencing. The value of 8% was empirically determined as tradeoff between detection of further variants and generating more false positive results. Further relaxing of the threshold did also not result in detection of all inspected variants due to stringent variant calling criteria ahead of this filtering step.

Further filtering of the variants in routine testing against minor allele frequency (MAF≤1%) and in silico predicted pathogenicity is performed as previously described [11] with nonsense, frameshift and canonical splice site variants being considered pathogenic.

For visualization of the identified SNVs, BAM-files were loaded into the SeqPilot SeqNext module (v4.1.2, JSI medical systems, Kippenheim, Germany) and IGV (Integrative Genomics Viewer) from Broad Institute.

All identified variants are validated by conventional Sanger sequencing. Divergent zygosity of variants may be corrected to the Sanger result after bioinformatic inspection of the alignment at the site of interest in IGV.

Copy number variation analysis

We performed copy number variation (CNV) analysis on highly covered samples sequenced on the Illumina Hiseq1500 system. Potential copy number alterations (CNA) were initially identified with the tools copynumber and copyCaller from VarScan v2.3.6 [33] on mapped reads with a maximum segment size of 300. All other parameter were used with standard settings. Thereby coverage of every target region of the sample of interest was internally normalized and compared versus normalized control data of other samples of the same run. CNVs were annotated using RefSeq gene file from UCSC (ftp://hgdownload.cse.ucsc.edu/ goldenPath/hg19/database/refGene.txt.gz). Potential CNVs were initially taken into account, if the CNV was detected by VarScan against at least 85% of the control patients and if the log₂ threshold was ≥ 0.6 in case of an amplification or ≤ −0.6 in case of a deletion. CNVs below this threshold (detected against 50% to 85% of the control patients) were indicated, but filtered against stringent criteria to be considered for subsequent MLPA validation: CNVs in critical exons with restricted discriminative mapping (see above and Fig. 2) and recurrent CNV artefacts in exons 43 and 46 were neglected. If heterozygous reads were evident in exons with an indicated CNV, these CNVs were also neglected as were CNVs in exons 1 and 42 due to insufficient or varying coverage. Results from CNV analysis were compared to previous results from MLPA (multiplex ligation dependent probe amplification) analysis. For the PKD1 gene the SALSA MLPA probemixes P351-C1 and P352-C2 were used (MRC-Holland, Amsterdam, The Netherlands). A total of six MLPA negative samples were used as negative controls for CNV detection.

As split read method we used the ClipCrop [34] program for identification of soft-clipped reads in PKD1. For patient 46, no CNV was explicitly called using the default thresholds of the program, but a manual review of the extracted split reads revealed an excess of such reads in a position in exon 15. These reads were visualized using IGV (Fig. 3). With 332 split reads supporting a breakpoint in that position, we further inspected the unmapped part of the split reads to find a potential matching deletion endpoint. A Smith-Waterman alignment of this sequence was performed versus the entire DNA sequence of PKD1, and, as suggested by coverage-based CNV analysis, the best match was found in the intronic region between exons 21 and 22. This represents a repeat-rich region that could not be enriched in the target capture process. Such a constellation might cause problems in the CNV calling module of ClipCrop as accompanying split reads at the potential deletion endpoint cannot be found. To further unravel the genomic structure of the patient at the CNV breakpoint, we applied the GATA tool [35] facilitating a more general motif search. We suppose that multiple re-shuffling of transposable repeat-elements in the genomically unstable regions of introns 21 and 22 could have created sequence alterations that cannot be mapped with conventional alignment programs. Application of the GATA tool decomposed the read sequence to short fragments in a sliding window approach (using windows of 20 bases) and subsequently mapped all fragments to the PKD1 sequence via BLASTN. The best-aligning fragments delimit the CNV region boundaries as shown in Fig. 3. A more detailed description of the used bioinformatic scripts are supplied in online S1 Materials And Methods.

Read simulation

In a first step, FASTA files for PKD1 exons 1–33 where extracted from the hg19 reference genome using exon coordinates from the RefSeq gene file (ftp://hgdownload.cse.ucsc.edu/goldenPath/hg19/database/refGene.txt.gz) as well as the FASTA files for PKD1 pseudogenes PKD1P1–6 using hg19 coordinates from GeneCards (e.g. http://www.genecards.org/cgi-bin/carddisp.pl?gene=PKD1P1) [36]. For all locations we added 20 bp of flanking sequence. Next, we calculated the length for every location and the amount of simulated reads to achieve an average coverage which will be equal to those of a MiSeq run under the assumption that we need 1000X coverage for every sequencing fragment with a length of 200 bp (S3 Fig.). After this, we simulated the 2 × 150 bp paired-end reads in FASTQ format with Wgsim from the SAMtools package (Li, H. wgsim—read simulator for next generation sequencing). For the parameter of Wgsim we used a base error rate of 0.001, an outer distance between the two ends of 200 bp and a length of the first and second read of 150 bp for all locations. For the reads of PKD1 exons 1–33 we simulated a mutation rate of 0.02 and for PKD1P1–6 a rate of 0 which simulates wild-type reads with a base error rate only and no mutations. The number of simulated paired-end reads was adjusted to the appropriate location. All other parameters were used with standard settings. After all reads had been simulated, we merged all simulated FASTQ files into one file for R1 and R2 respectively and analyzed them with our in-house bioinformatic pipeline (see Bioinformatic data analysis). For more details see S1 Materials And Methods.

Ultra-deep sequencing

For only targeting the PKD1 coding region, long-range PCR (LR-PCR) amplicons were generated as described elsewhere using locus-specific primers [8,9]. Concentration of amplicons was measured by a fluorescence-based method. To minimize variation of coverage between amplified exons, LR-amplicons were equimolarly pooled or based on empirically determined ratios. These amplicon pools were taken as input for library preparation using the NEBNext DNA library prep master mix set for Illumina (New England Biolabs, Ipswich, MA, USA) as described above. Libraries of five to six patients were pooled and subjected to sequencing on an Illumina MiSeq system (2 × 150 bp PE reads). Data analysis was performed by the standard bioinformatics pipeline.

Data repository

Sequence data for PKD1 have been deposited at the European Genome-phenome Archive (EGA), which is hosted by the EBI, under accession number EGAS00001001003. Access can be gained upon request to the authors.