doi: 10.1111/ahg.12060

Effective Variant Detection by Targeted Deep Sequencing of DNA Pools: An Example from Parkinson’s Disease Lasse Pihlstrøm1 ∗ , Aina Rengmark1 , Kari Anne Bjørnara˚ 2 and Mathias Toft1 1 Oslo University Hospital, Department of Neurology, P.O. Box 4950 Nydalen, 0424 Oslo, Norway 2 Drammen Hospital, Department of Neurology, 3004 Drammen, Norway

Summary Next-generation sequencing technologies will dominate the next phase of discoveries in human genetics, but considerable costs may still represent a limitation for studies involving large sample sets. Targeted capture of genomic regions may be combined with deep sequencing of DNA pools to efficiently screen sample cohorts for disease-relevant mutations. We designed a 200 kb HaloPlex kit for PCR-based capture of all coding exons in 71 genes relevant to Parkinson’s disease and other neurodegenerative disorders. DNA from 387 patients with Parkinson’s disease was combined into 39 pools, each representing 10 individuals, before library preparation with barcoding and Illumina sequencing. In this study, we focused the analysis on six genes implicated in Mendelian Parkinson’s disease, emphasizing quality metrics and evaluation of the method, including validation of variants against individual genotyping and Sanger sequencing. Our data showed 97% sensitivity to detect a single nonreference allele in pools, rising to 100% where pools achieved sequence depth above 80x for the relevant position. Pooled sequencing detected 18 rare nonsynonymous variants, of which 17 were validated by independent methods, corresponding to a specificity of 94%. We argue that this design represents an effective and reliable approach with possible applications for both complex and Mendelian genetics. Keywords: Next-generation sequencing, pooled sequencing, targeted NGS, mutation screening, Parkinson’s disease

Introduction An increasingly available and affordable technology, nextgeneration sequencing (NGS) is now becoming the primary discovery tool in human genetics (Goldstein et al., 2013). This development also calls for rational, cost-effective and scalable study designs to generate sequence data for different scientific applications. Several early studies have successfully used an NGS approach to identify mutations underlying rare Mendelian disorders (Lupski et al., 2010; Ng et al., 2010), and many geneticists are now looking towards larger sample sets and common phenotypes. Sequence data will be crucial to identify the functionally relevant variation behind association signals from genomewide association studies (GWAS) and further explore the genetic architecture of common diseases. For genetically het∗

Corresponding author: Lasse Pihlstrøm, Department of Neurology, Oslo University Hospital, P.O. Box 4950 Nydalen, N-0424 Oslo, Norway. Tel: +47 23079023; Fax: +47 23070490; E-mail: [email protected]

 C

2014 John Wiley & Sons Ltd/University College London

erogeneous conditions, NGS might also be the most rational approach to screen for known mutations. Individual wholegenome or whole-exome sequencing is currently being used to study complex genetics, but will in many cases still be too resource-demanding for large sample sets. In this article, we outline a study design combining targeted capture with deep sequencing of DNA pools and demonstrate how this strategy may be applied to effectively sequence relevant genes in 387 patients with Parkinson’s disease (PD). PD is a common neurodegenerative disorder of complex aetiology. To date, 20 genetic loci influencing PD susceptibility have been identified through GWAS (Nalls et al., 2011; Lill et al., 2012; Pankratz et al., 2012). However, a minority of patients have monogenic forms of the disease, following dominant or recessive inheritance patterns. The proportion of PD caused by mutations in Mendelian genes varies considerably between populations (Nuytemans et al., 2010; Puschmann, 2013). Where previous efforts to determine mutation frequencies in large sample cohorts have been limited to a few mutations or genes, NGS now provides the opportunity for large-scale screening strategies.

Annals of Human Genetics (2014) 78,243–252

243

L. Pihlstrøm et al.

We designed a targeted DNA capture panel of 997 exons in 71 genes relevant to PD. These included all proposed and established Mendelian genes, 24 genes adjacent to GWAS top hits, as well as a range of genes implicated in neurodegenerative disorders with overlapping clinical or pathological features. This panel was used for pooled sequencing of DNA from 387 PD patients in order to generate a data set that may serve as a resource for several different subprojects. The strategy for downstream filtering and validation of variants will need to be attuned to the specific hypothesis. In this study, our primary aim was to evaluate the performance of the targeted pooled NGS design. For this purpose, we decided to focus on the subset of genes where the coding variability in PD is currently best characterized, namely those implicated in Mendelian disease. We defined a target region corresponding to all coding regions of the six well-established autosomal dominant (SNCA; MIM# 163890, LRRK2; MIM# 609007 and VPS35; MIM# 601501) and recessive (PARK2; MIM# 602544, PINK1; MIM# 608309 and PARK7 or DJ-1; MIM# 602533) PD genes. In this article, we present the results from analysis of these genes with emphasis on quality measures and methodological issues.

Materials and Methods Subjects The study was approved by the Regional Committee for Medical Research Ethics (Oslo, Norway). All participants gave written, informed consent. Patients fulfilling United Kingdom Parkinson’s Disease Brain Bank criteria for PD were recruited between 2007 and 2012 at Oslo University Hospital and Drammen Hospital, Norway. The mean age of onset was 53 years (range 27–78) and the proportion of males was 66%. All participants were examined by a neurologist and blood samples were collected by venopuncture.

Preparation of DNA Pools DNA was extracted from whole blood by the same standard techniques for all samples. A selection of samples was tested on 1% agarose gel to assure the integrity of DNA. DNA was homogenized for 30 min in a thermoshaker at 50◦ C, and all samples were diluted to a working solution of approximately 50 ng/μl. Each sample was then carefully measured on a Qubit 2.0 Fluorometer (Invitrogen/Life Technologies, Paisley, UK) and further diluted with TE-buffer to 20 ng/μl. Subsequently, 10 μl, corresponding to 200 ng of DNA, were drawn from each of the samples, and mixed together with others in pools representing 10 individuals. We created 39

244

Annals of Human Genetics (2014) 78,243–252

pools from a total of 387 individuals, three samples being present in two different pools.

Capture, Enrichment and Barcoding The HaloPlex Target Enrichment System (Agilent Technologies, Santa Clara, CA, USA) relies on a tailored cocktail of restriction enzymes and customized probes to capture genomic regions of interest, which are subsequently amplified by polymerase chain reaction (PCR). We used the HaloPlex online design tool to create a 200 kb panel targeting 997 exons comprising 71 genes relevant to PD (Customized HaloPlex Kit 48 rxn, Art no.96005, Halo Genomics, currently Agilent Technologies, https://earray.chem.agilent.com/suredesign/). We used 39 unique barcodes from the HaloPlex Kit for this experiment. Target enrichment was performed according to the manufacturer’s protocol (HaloPlex Target Enrichment System for Illumina Sequencing Version A, February 2012 and HaloPlex PCR Target Enrichment & Library Preparation Guide version 2.0, November 2011). In brief, 45 μl aliquots of each pool, corresponding to 900 ng of genomic DNA, were digested by restriction enzymes. Successful digestion was verified by gel electrophoresis, demonstrating expected bands at 125, 175 and 475 bp. DNA pools were then hybridized overnight to customized, biotinylated probes. Hybridized probes were captured with magnetic beads and target fragments were ligated to create circular DNA molecules. Subsequently, libraries were amplified by PCR, introducing unique index sequences that allow all pools to be sequenced together.

Sequencing All libraries of target-enriched pooled DNA were analysed on a Bioanalyzer 2100 with high-sensitivity DNA chips (Agilent Technologies) to verify successful enrichment, demonstrating a smear of amplicons ranging from 225 to 525 bp. The samples were sequenced at the Norwegian Sequencing Centre, Oslo. All samples were run together on two lanes of an Illumina HiSeq 2000 to perform 100 bp paired-end sequencing.

Bioinformatic Analysis of Sequencing Data A detailed description of the bioinformatic pipeline indicating the specific tools and parameters used is provided in Table S1. The HaloPlex workflow may create target fragments that are short enough for 100 bp sequencing to read into the adapter towards the end of reads. Adapter sequence and low-quality tails of reads were therefore removed with the software Trimmomatic 0.30 (Lohse et al., 2012). Reads were aligned to the

 C

2014 John Wiley & Sons Ltd/University College London

Targeted Pooled Sequencing in PD

reference genome (GRCh37) with bwa 0.5.9 (Li & Durbin, 2009). Picard 1.77 was used to index and compress aligned sequence files to bam-format. Further sequence data processing, assessment of coverage and mismatch rates, variant calling and filtration was performed with the Genome Analysis Toolkit, GATK 2.5.2 (DePristo et al., 2011), restricting the target region to PARK7 (NM_007262.4), PINK1 (NM_032409.2), SNCA (NM_000345.3), PARK2 (NM_004562.2), LRRK2 (NM_198578.3) and VPS35 (NM_018206.4). We used ANNOVAR (version 2012may25) to annotate called variants (Wang et al., 2010). The error rate of NGS usually requires some kind of quality filtering to avoid large numbers of false-positive variants. For this study, we applied quality filters with loose parameters in order to maintain a high sensitivity (Table S1). Restricting the analysis to mutations of potential relevance to monogenic PD, we further filtered out synonymous single-nucleotide polymorphisms (SNPs) and mutations that are classified as “not pathogenic” in the Parkinson Disease Mutation Database (Cruts et al., 2012). In a final filtering step, we removed variants with a frequency above 0.03 in the Exome Variant Server (ESP4500) or 1000 genomes (1000g2012feb) databases (Fig. 1).

Positive Controls and Validation To serve as positive controls, 11 exonic variants in SNCA, PARK2, PINK1 and LRRK2 were genotyped prior to sequencing in 350 individual samples by MALDI-TOF mass spectrometry using the Sequenom MassARRAY system. For all variants passing filters, we validated mutations and identified individual carriers by Sanger sequencing of all samples represented in the positive pools. We used Primer3Plus to design primers for PCR amplification and sequencing of relevant exons. Primers and conditions are available on request. Amplicons were sequenced bidirectionally and Sequencher 5.1 (Gene Codes Corporation, Ann Arbor, MI, USA) software was used for data analysis.

Results Coverage and Mismatch Rate The six monogenic PD genes we analysed in this study comprise a total of 99 exons, corresponding to a target region of 14 kb. The total coverage across this region was 2.5 Gb, which equals each position on average being sequenced to a depth of approximately 4600x for each of the pools. Benchmarks used when assessing coverage in pooled sequencing vary considerably between previously published studies. As described below, we did not observe missed alleles when depth was above

 C

2014 John Wiley & Sons Ltd/University College London

Figure 1 Overview of variant filtering. The figure indicates how 60 initially called variants were processed through a series of filtering steps to a final list of 17 validated, potentially pathogenic variants.

80x, corresponding to each allele being read on average four times. The proportion of bases covered above 80x in 80% of pools was 95% for the whole target region and above 94% for each individual gene. The performance of each pool is illustrated in Figure 2. Setting a stricter coverage benchmark, at least 200x or 10 reads of each allele on average, decreases these measures only slightly, to 93% of the total target and above 92% for each gene. We note that one out of 39 pools achieved markedly less coverage than the rest of the samples, with only 79% of the target region covered at 80x (Fig. 2, Pool 18). The same pool was also characterized as having the lowest library quality as analysed by Bioanalyzer prior to sequencing. Next, we assessed the occurrence of complete coverage gaps in the data. The best performing pools had reads mapping to 98.8% of the target region, corresponding exactly to the theoretically predicted coverage of the HaloPlex kit (data not shown). Coverage gaps occurred mainly in LRRK2 exons 9, 20 and 32, and in VPS35 exon 14 (Table 1).

Annals of Human Genetics (2014) 78,243–252

245

L. Pihlstrøm et al.

Figure 2 Coverage of target region across pools. The figure shows the proportion of targeted exonic positions covered above 80x depth for each of the 39 pools.

Pooled sequencing relies on the ability to distinguish a true allele occurring on a small proportion of reads from the background error rate of the experiment. In Illumina 100 bp sequencing, the quality tends to fall towards the end of reads. This might make a subset of targeted positions vulnerable to higher error rates, given that the PCR-based HaloPlex protocol produces many identical reads, with only few overlaps per position. We therefore applied strict trimming of low-quality ends of reads in the raw data (Table S1). Furthermore, we assessed the mismatch rate, meaning the proportion of bases differing from reference, calculated across each instrument cycle. There was no systematic relationship between mismatch

rate and read position (data not shown). The overall mismatch rate for the experiment was 0.13%.

Sensitivity To assess the sensitivity of the pooled sequencing experiment, we compared results from genotyping of 350 patients with called variants from sequencing of the corresponding 35 pools. A detailed account of this comparison is given in Table S2. We genotyped 11 exonic SNPs, but for three of these, all patients were homozygous for the reference allele. For the

Table 1 Depth and coverage statistics. Gene

PARK7

PINK1

SNCA

PARK2

LRRK2

VPS35

Total

Percentage of exonic positions theoretically targeted by HaloPlex capture design Percentage of total exonic positions sequenced to minimum 80x in 80% of samples Average depth per exonic position per sample pool Number of exons Exons flagged for coverage gaps1

100

99.8

100

99.5

98.5

98.2

98.8

100

97.7

100

95.1

95.6

94.8

95.2

5627 6

4176 8

5889 5

4581 12

4928 51 Exons 9, 20 and 32

3570 17 Exon 14

4627 99

1

246

The GATK tool DiagnoseTargets was used to flag exons if more than 20% of pools lacked coverage for 10% of positions.

Annals of Human Genetics (2014) 78,243–252

 C

2014 John Wiley & Sons Ltd/University College London

Targeted Pooled Sequencing in PD

remaining eight variants, we observed heterozygous carriers, with the total number of nonreference alleles in the data set ranging from 1 to 205 for each SNP. The most fundamental aspect of sensitivity in this design will be the experiment’s ability to detect a variant where only one allele is present in a pool. We tested 59 such single-allele pools representing seven different SNPs and a variant was called in 57 cases, corresponding to a sensitivity of 97%. The two cases where sequencing failed concerned rs55774500 in exon 3 of PARK2, and depth was below 80x at this position for both pools where the variant was missed. Consequently, if we restrict the analysis to pools achieving a depth above 80x at the relevant position, the observed sensitivity was 100%. To validate findings, we performed a total of 210 individual bidirectional Sanger sequencing analyses across 12 different exons. We observed no additional variants missed by pooled NGS in this Sanger sequencing data, amounting to a total of more than 37 kb. Where a pool contains multiple nonreference alleles for a given SNP, the variant calling algorithm will calculate the most probable allele count. For studies concerned with rare variant detection, the accuracy of this estimate is less important, especially if positive pools are reanalysed by Sanger sequencing. However, in study designs where pooled sequencing is used to estimate allele frequencies of common variants, this quantitative aspect is crucial. We tested 56 instances where multiple nonreference alleles of an SNP, from 2 up to 10, were present in a pool. We observed that the number of nonreference alleles was called correctly in a majority of cases, but could vary with up to two alleles more or less as compared to genotyping (Table S2). The precision tended to fall with a higher proportion of nonreference alleles in a tested pool. No variants were called where pools were negative for an SNP as assessed by genotyping.

Filtering and Detection of Variants Initial variant calling identified 60 variants in the data set passing the basic quality threshold. Subsequently, we performed a series of filtering steps retaining 18 nonsynonymous, low-frequency variants of potential pathogenic significance (Fig. 1). Out of these 18 variants, two (PARK2 p.Ala82Glu and LRRK2 p.Gly2019Ser) were already genotyped as positive controls, and one (LRRK2 p.Asn1437His) was found to have been previously detected in a diagnostic setting. Sanger sequencing confirmed 14 out of the 15 remaining variants and individual carriers were identified in each of the implicated pools. One single variant was not validated by Sanger sequencing and discarded as a false positive, corresponding to an observed specificity of 94% in our experiment. The final set of 17 variants is listed in Table 2. Based on presence in db-

 C

2014 John Wiley & Sons Ltd/University College London

SNP137, Exome Variant Server and 1000 genomes databases, 11 variants were previously reported and 6 were novel. The novel variants were submitted to the ClinVar database (Landrum et al., 2014). We note that although we chose to filter out synonymous SNPs in this study, this class of variants may also contribute to disease through mechanisms such as splicing and RNA processing (Sauna & Kimchi-Sarfaty, 2011). Such effects are difficult to predict for variants detected in a screening context, but would have to be considered in a complete analysis of possibly pathogenic mutations.

Discussion In this study, we performed pooled deep sequencing of 387 PD samples, using HaloPlex PCR-based target enrichment to capture relevant genomic regions. We identified 11 known and 6 novel variants in six Mendelian PD genes and evaluated the method by coverage assessment and validation against SNP genotyping and Sanger sequencing. We suggest that the combination of targeted capture and DNA pooling may be a cost-effective and scalable approach for variant detection in large sample sets in a research setting. In Table 3, we summarize key aspects of a pooled NGS design as compared to a traditional Sanger sequencing approach. A targeted pooled sequencing strategy was recently used successfully to identify rare causal variants in breast and ovarian cancer (Ruark et al., 2013), as well as short stature (Wang et al., 2013). These studies investigated large numbers of candidate genes, more than 1000 in the latter case. In contrast, a study of Alzheimer’s disease applied pooled sequencing to screen for mutations in only five genes (Cruchaga et al., 2012). We designed a panel sufficiently large to entail all genes implicated in previous genetic studies of PD and related disorders, but not extensive “new” candidate gene sets. The main motivation for creating DNA pools is to reduce the costs related to target enrichment and sequencing, which still represent a substantial barrier for larger projects. While some laboratory effort is required for the preparation of equimolar pools, the manual workload is considerably reduced downstream, as compared to individual sequencing. A major advantage of the HaloPlex-targeted enrichment system is the convenient laboratory workflow, integrating both capture and library preparation. The protocol allows one person to prepare a set of finished libraries within two working days, and requires no larger specialized laboratory instruments. Barcoding offers the opportunity to mix DNA libraries on the same flow cell lane, while allowing raw data to be extracted for each individual library. We note, however, that our pooled design maximizes efficiency on two different levels by sequencing barcoded pools, which reduces the experiment

Annals of Human Genetics (2014) 78,243–252

247

L. Pihlstrøm et al.

Table 2 Results from targeted pooled sequencing of Mendelian PD genes. Genomic coordinate

Gene1

Base change

Amino acid change

dbSNP137

ClinVar submission

1:8037783 1:8037788 1:20964591 1:20971129 1:20971158 1:20975105 6:161771219 6:162394367 6:162683724 12:40631852 12:40677787 12:40688681 12:40692281 12:40703027 12:40734202 12:40758762 16:46702913

PARK7 PARK7 PINK1 PINK1 PINK1 PINK1 PARK2 PARK2 PARK2 LRRK2 LRRK2 LRRK2 LRRK2 LRRK2 LRRK2 LRRK2 VPS35

c.394A>G c.399G>C c.644C>T c.923T>A c.952A>T c.1231G>A c.1310C>T c.701G>A c.245C>A c.518A>G c.2352C>A c.2843G>A c.3333G>T c.4309A>C c.6055G>A c.7300A>G c.1576C>T

p.Lys132Glu p.Met133Ile p.Pro215Leu p.Leu308Gln p.Met318Leu p.Gly411Ser p.Pro437Leu p.Arg234Gln p.Ala82Glu p.Asn173Ser p.Ser784Arg p.Arg948Gln p.Gln1111His p.Asn1437His p.Gly2019Ser p.Ile2434Val p.Arg526Cys

rs200894731

ESP database frequency2

SCV000114936 SCV000114938 SCV000114939 rs139226733 rs45478900 rs149953814 rs144032774 rs55774500 rs201415714

0.00093 0.001022 0.002138 0.000186 0.00158 SCV000114940

rs200442352 rs78365431 rs74163686 rs34637584

0.000093 0.000465 SCV000114941 SCV000114937

Allele count3 1 1 1 1 3 1 4 1 16 1 1 1 1 1 3 1 1

The table lists our final set of validated variants of potential pathogenic significance detected in 387 patients by pooled resequencing of Mendelian PD genes. 1 Genomic coordinates refer to build GRCh37. Relevant RefSeq accession numbers for mRNA and protein are: PARK7: NM_007262.4, NP_009193.2; PINK1: NM_032409.2, NP_115785.1; PARK2: NM_004562.2, NP_004553.2; LRRK2: NM_198578.3, NP_940980.3; VPS35: NM_018206.4, NP_060676.2. 2 ESP database frequencies are given as proportions of 1. 3 All mutation carriers were heterozygotes. Table 3 Comparison of study designs for rare variant detection in disease cohorts. Individual Sanger sequencing

Pooled targeted next-generation sequencing

Project may develop over time Requires a series of experiments to analyse individual exons in individual samples Accumulating costs and workload with increasing data generation Large-scale projects very labour-intensive and time-consuming

Highly accurate results, method also applicable in a diagnostic setting Data analysis straightforward with conventional tools Validation of findings generally not required

Project relies on the initial design Exons and samples in the range of hundreds, or even thousands, are analysed in parallel in a single experiment One initial investment effectively generates large amounts of data Easily scaled to projects involving large sample sets and extensive target regions Target enrichment may be incomplete in challenging genomic regions Sensitivity and specificity adequate in the context of clinical research Requires more advanced bioinformatic processing of data May require validation of findings by independent methods

costs considerably more than barcoding alone. Savings are most significant for the target enrichment, which is currently the most expensive step in this kind of study, but efficiency is also further optimized for the sequencing, utilizing the whole capacity of the instrument run. For our 200 kb target, ample coverage for all samples was obtained on only two lanes of the Illumina HiSeq.

We achieved coverage above 80x for 96% of exonic positions in 80% of samples, compared to 98.8% theoretically targeted by our capture kit. This corresponds to high performance regarding both design efficiency and target enrichment efficiency, as compared to previous successful pooled sequencing studies (Rivas et al., 2011; Ruark et al., 2013). By comparison, frequently used capture technologies based on random

PCR primers successively tailored to amplify full target region

248

Annals of Human Genetics (2014) 78,243–252

 C

2014 John Wiley & Sons Ltd/University College London

Targeted Pooled Sequencing in PD

DNA shearing and hybridization will typically miss considerably more of the region of interest due to limitations in design, even though deep sequencing yields good coverage of the actual baits in the kit (Hedges et al., 2011; Kiialainen et al., 2011). Bearing this in mind, we also recognize that the target capture in our study has limitations. HaloPlex technology relies on restriction sites and customized hybridizing probes, which may be difficult to tailor for some genomic regions, thus limiting design efficiency. Our data show that while some reads were present for the total theoretical target, a minor proportion of the region failed to achieve sufficient coverage. Considering that sequence depth per position in each pool was 4600x on average, and that varying the coverage benchmark made little impact on the statistics, it is unlikely that increasing the total sequencing depth would rescue these poorly covered fragments. The phenomenon of variable relative coverage is well described in the NGS literature (Ross et al., 2013). Loci with extreme base compositions (i.e., highly GC- or AT-rich regions) are known to be the major cause of low relative coverage, where library amplification by PCR has been shown to be the most critical step (Aird et al., 2011). In our data, we observed that a small subset of pools, and one in particular, showed markedly lower coverage performance than the majority of pools. This suggests that factors such as subtle variation in the quality of input DNA or the inevitable methodological variability introduced by manual processing of samples may impact the enrichment efficiency in these most vulnerable target regions. Validation against positive control genotypes demonstrated how stretches of inadequate coverage might lead to false negatives, as our pooled NGS experiment missed two alleles where sequence depth was low. However, if the analysis is restricted to adequately covered loci, our observed sensitivity was 100%. This indicates that the principle of pooling DNA from 10 individuals prior to target enrichment is compatible with excellent sensitivity, whereas vulnerability to false negatives primarily depends on the enrichment efficiency. While this study is concerned with effective mutation screening in large sample sets for research purposes, the required level of sensitivity would be higher in the context of clinical diagnostics. A recent study reported the performance of individual targeted NGS for clinical diagnosis of hereditary cardiomyopathy (Sikkema-Raddatz et al., 2013). Achieving 99% coverage in their experiment, the authors recommend that diagnostic NGS should be supplemented by Sanger sequencing where relevant genes are inadequately covered. A similar approach might also be considered for mutation screening in a research context, if sensitivity is to be increased further.

 C

2014 John Wiley & Sons Ltd/University College London

Regarding the specificity, the pooled design relies on accurate distinction of true nonreference alleles occurring in a low proportion of reads from the background of sequencing errors. The empirical mismatch rate entails all deviations from reference, both biological variation and miscalls of technical origin, although sequencing errors will dominate quantitatively in highly conserved regions where true variants are rare. We observed a total mismatch rate of 0.13% in our experiment. This figure corresponds well to reported estimations of error rate in Illumina sequencing (Glenn, 2011), and is far below the expected 5% read frequency of a single allele in our pools of 10 individuals. It has been shown, however, that the distribution of errors is not random, but tends to follow specific patterns (Nakamura et al., 2011). While such sequence-specific susceptibility to errors may give rise to false-positive variant calls, the data quality metrics will typically be low in such loci, allowing for bioinformatic filtering to increase specificity. Much like for NGS in general, there is currently no standardized pipeline or established consensus on the optimal bioinformatic strategies for the analysis of pooled sequencing data. It may be worth noting that restriction enzyme digestion and PCR produces sequence data with many identical reads, in contrast to reads from randomly fragmented DNA. This will affect some of the quality-control metrics that are often exploited in variant filtering. The parameters and strategies for quality filtering will have to be tuned to the relevant hypothesis and design. For this study, we wanted to prioritize a high sensitivity and we therefore chose loose parameters that filtered out only a few variants. Still, all but one variant were confirmed by Sanger sequencing, corresponding to a specificity of 94%. For other pooled sequencing applications, one might be more concerned with high specificity and apply stricter filtering. A weakness in our study design is the need to follow up positive findings by Sanger sequencing of all 10 individuals in a pool. We note, however, that this step may not be necessary for all pooled NGS experiments. As an example, a recent study of rheumatoid arthritis used pooled sequencing to identify rare risk variants in a design that relied on high-quality calls without validation or identification of individual carriers (Diogo et al., 2013). The Genome Analysis Toolkit provides a comprehensive framework for sequence data processing and variant detection that will be familiar to many researchers involved in NGS analysis. Since the GATK version 2.0, an extension to the UnifiedGenotyper tool has been implemented to support variant calling from pooled samples. Noting that it represents only one among a range of available software tools (Bansal, 2010; Rivas et al., 2011; Vallania et al., 2012), we would argue that our results demonstrate the capacity of UnifiedGenotyper to effectively and reliably call SNPs from pooled sequencing.

Annals of Human Genetics (2014) 78,243–252

249

L. Pihlstrøm et al.

A clear weakness concerning this study is that our set of positive control genotypes did not include any small indels. Indel calling in NGS is generally more challenging than SNP calling, and could be expected to be even more difficult in the context of pooled sequencing. The detection of larger copy number variations (CNVs) would require different methods. We note that for PARK2 and SNCA in particular, a substantial proportion of pathogenic mutations are CNVs (Nuytemans et al., 2010). With the exception of possible undetected CNVs, our experiment gives an overview of potentially pathogenic Mendelian PD gene mutations in a sample set of 387 patients. While we identified 17 variants from 38 different carriers, only a small proportion of these will have monogenic PD. For most of the listed variants, the pathogenic potential is currently uncertain. Twenty-eight patients had heterozygous mutations in recessive genes PARK7, PINK1, or PARK2. It remains unclear to what extent the carrier state of such mutations represents a risk factor for PD (Nuytemans et al., 2010). Regarding the dominant genes, four patients had highly penetrant LRRK2 mutations (three p.Gly2019Ser, one p.Asn1437His) (Ross et al., 2011). The remaining LRRK2 variants have an unknown role in relation to PD. Coding variants in SNCA are very rare, but recently a few novel pathogenic mutations have been reported (Appel-Cresswell et al., 2013; Lesage et al., 2013; Proukakis et al., 2013). We found no SNCA variants in our sample set. From our final set of 17 validated variants, 6 were novel mutations not previously reported. This demonstrates the capacity of the pooled sequencing approach to identify new variants. One novel variant was found in VPS35. This gene was implicated in dominant PD more recently (Vilarino-Guell et al., 2011; Zimprich et al., 2011). The p.Arg526Cys mutation is classified as damaging by prediction tools SIFT and MutationTaster, but not by PolyPhen2, and was detected in a patient without known family history of PD. It remains to be seen whether this variant is found in other individuals with or without PD. In this study, we analysed sequence data from 6 genes included in a targeted enrichment panel of 71 genes in total. Having gained experience with the method and demonstrated its ability to reliably detect rare coding variants, we will make use of the remaining sequence data in future projects to explore various hypotheses related to the genetic architecture of Parkinsonism. In summary, we have demonstrated how deep sequencing of genomic target regions in DNA pools may be a costeffective, scalable and reliable design for NGS studies in a research context. The example of PD, a complex disorder where a minority of patients have Mendelian forms of the disease, serves to show how targeted pooled sequencing can be

250

Annals of Human Genetics (2014) 78,243–252

useful to pursue a range of scientific hypotheses, concerning both complex and monogenic conditions.

Acknowledgements L. Pihlstrøm is supported by a grant from South-Eastern Norway Regional Health Authority, Norway. M. Toft and A. Rengmark are supported by grants from the Research Council of Norway. Patient recruitment and sample collection has also been funded by the Norwegian Parkinson Research Fund and Reberg’s Legacy. The authors declare no conflicts of interest. The authors thank research nurse Lena Pedersen, Department of Neurology, Oslo University Hospital for assistance in sample collection and handling, Zafar Iqbal, Department of Neurology, Oslo University Hospital for comments to the manuscript, and Tim Hughes, the Norwegian Sequencing Centre, for bioinformatic advice. The sequencing service was provided by the Norwegian Sequencing Centre (http://www.sequencing.uio.no), a national technology platform hosted by the University of Oslo and supported by the “Functional Genomics” and “Infrastructure” programs of the Research Council of Norway and the South-Eastern Norway Regional Health Authority.

Web Resources The URLs for databases and software referenced in this article are as follows: Online Mendelian Inheritance in Mam (OMIM), http:// omim.org 1000 Genomes, http://www.1000genomes.org NHLBI Exome Sequencing Project (ESP), http://evs.gs. washington.edu/EVS/ ClinVar, http://www.ncbi.nlm.nih.gov/clinvar/ dbSNP, http://www.ncbi.nlm.nih.gov/projects/SNP/ Burrows-Wheeler Aligner, http://bio-bwa.sourceforge. net Picard, http://picard.sourceforge.net The Genome Analysis Toolkit, http://www.broadinstitute. org/gatk/ Parkinson Disease Mutation Database, http://www. molgen.ua.ac.be/PDmutDB/ ANNOVAR, http://www.openbioinformatics.org/ annovar/ Trimmomatic, http://www.usadellab.org/cms/?page= trimmomatic SIFT, http://sift.jcvi.org PolyPhen2, http://genetics.bwh.harvard.edu/pph2/ MutationTaster, http://www.mutationtaster.org Primer3Plus, http://primer3plus.com

 C

2014 John Wiley & Sons Ltd/University College London

Targeted Pooled Sequencing in PD

References Aird, D., Ross, M. G., Chen, W. S., Danielsson, M., Fennell, T., Russ, C., Jaffe, D. B., Nusbaum, C. & Gnirke, A. (2011) Analyzing and minimizing PCR amplification bias in Illumina sequencing libraries. Genome Biol 12, R18. Appel-Cresswell, S., Vilarino-Guell, C., Encarnacion, M., Sherman, H., Yu, I., Shah, B., Weir, D., Thompson, C., Szu-Tu, C., Trinh, J., Aasly, J. O., Rajput, A., Rajput, A. H., Jon Stoessl, A. & Farrer, M. J. (2013) Alpha-synuclein p.H50Q, a novel pathogenic mutation for Parkinson’s disease. Mov Disord 28, 811–813. Bansal, V. (2010) A statistical method for the detection of variants from next-generation resequencing of DNA pools. Bioinformatics 26, i318–i324. Cruchaga, C., Haller, G., Chakraverty, S., Mayo, K., Vallania, F. L., Mitra, R. D., Faber, K., Williamson, J., Bird, T., Diaz-Arrastia, R., Foroud, T. M., Boeve, B. F., Graff-Radford, N. R., St Jean, P., Lawson, M., Ehm, M. G., Mayeux, R. & Goate, A. M. (2012) Rare variants in APP, PSEN1 and PSEN2 increase risk for AD in late-onset Alzheimer’s disease families. PLoS One 7, e31039. Cruts, M., Theuns, J. & Van Broeckhoven, C. (2012) Locus-specific mutation databases for neurodegenerative brain diseases. Hum Mutat 33, 1340–1344. DePristo, M. A., Banks, E., Poplin, R., Garimella, K. V., Maguire, J. R., Hartl, C., Philippakis, A. A., Del Angel, G., Rivas, M. A., Hanna, M., Mckenna, A., Fennell, T. J., Kernytsky, A. M., Sivachenko, A. Y., Cibulskis, K., Gabriel, S. B., Altshuler, D. & Daly, M. J. (2011) A framework for variation discovery and genotyping using next-generation DNA sequencing data. Nat Genet 43, 491–498. Diogo, D., Kurreeman, F., Stahl, E. A., Liao, K. P., Gupta, N., Greenberg, J. D., Rivas, M. A., Hickey, B., Flannick, J., Thomson, B., Guiducci, C., Ripke, S., Adzhubey, I., Barton, A., Kremer, J. M., Alfredsson, L., Sunyaev, S., Martin, J., Zhernakova, A., Bowes, J., Eyre, S., Siminovitch, K. A., Gregersen, P. K., Worthington, J., Klareskog, L., Padyukov, L., Raychaudhuri, S. & Plenge, R. M. (2013) Rare, low-frequency, and common variants in the protein-coding sequence of biological candidate genes from GWASs contribute to risk of rheumatoid arthritis. Am J Hum Genet 92, 15–27. Glenn, T.C. (2011) Field guide to next-generation DNA sequencers. Mol Ecol Resour 11, 759–769. Goldstein, D. B., Allen, A., Keebler, J., Margulies, E. H., Petrou, S., Petrovski, S. & Sunyaev, S. (2013) Sequencing studies in human genetics: design and interpretation. Nat Rev Genet 14, 460–470. Hedges, D. J., Guettouche, T., Yang, S., Bademci, G., Diaz, A., Andersen, A., Hulme, W. F., Linker, S., Mehta, A., Edwards, Y. J., Beecham, G. W., Martin, E. R., Pericak-Vance, M. A., Zuchner, S., Vance, J. M. & Gilbert, J. R. (2011) Comparison of three targeted enrichment strategies on the SOLiD sequencing platform. PLoS One 6, e18595. Kiialainen, A., Karlberg, O., Ahlford, A., Sigurdsson, S., LindbladToh, K. & Syvanen, A. C. (2011) Performance of microarray and liquid based capture methods for target enrichment for massively parallel sequencing and SNP discovery. PLoS One 6, e16486. Landrum, M. J., Lee, J. M., Riley, G. R., Jang, W., Rubinstein, W. S., Church, D. M. & Maglott, D. R. (2014) ClinVar: public archive of relationships among sequence variation and human phenotype. Nucleic Acids Res 42, D980–D985. Lesage, S., Anheim, M., Letournel, F., Bousset, L., Honore, A., Rozas, N., Pieri, L., Madiona, K., Durr, A., Melki, R., Verny, C. & Brice, A. (2013) G51D alpha-synuclein mutation causes

 C

2014 John Wiley & Sons Ltd/University College London

a novel parkinsonian-pyramidal syndrome. Ann Neurol 73, 459– 471. Li, H. & Durbin, R. (2009) Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics 25, 1754–1760. Lill, C. M., Roehr, J. T., Mcqueen, M. B., Kavvoura, F. K., Bagade, S., Schjeide, B. M., Schjeide, L. M., Meissner, E., Zauft, U., Allen, N. C., Liu, T., Schilling, M., Anderson, K. J., Beecham, G., Berg, D., Biernacka, J. M., Brice, A., Destefano, A. L., Do, C. B., Eriksson, N., Factor, S. A., Farrer, M. J., Foroud, T., Gasser, T., Hamza, T., Hardy, J. A., Heutink, P., Hill-Burns, E. M., Klein, C., Latourelle, J. C., Maraganore, D. M., Martin, E. R., Martinez, M., Myers, R. H., Nalls, M. A., Pankratz, N., Payami, H., Satake, W., Scott, W. K., Sharma, M., Singleton, A. B., Stefansson, K., Toda, T., Tung, J. Y., Vance, J., Wood, N. W., Zabetian, C. P., Young, P., Tanzi, R. E., Khoury, M. J., Zipp, F., Lehrach, H., Ioannidis, J. P. & Bertram, L. (2012) Comprehensive research synopsis and systematic meta-analyses in Parkinson’s disease genetics: The PDGene database. PLoS Genet 8, e1002548. Lohse, M., Bolger, A. M., Nagel, A., Fernie, A. R., Lunn, J. E., Stitt, M. & Usadel, B. (2012) RobiNA: a user-friendly, integrated software solution for RNA-Seq-based transcriptomics. Nucleic Acids Res 40, W622–W627. Lupski, J. R., Reid, J. G., Gonzaga-Jauregui, C., Rio Deiros, D., Chen, D. C., Nazareth, L., Bainbridge, M., Dinh, H., Jing, C., Wheeler, D. A., Mcguire, A. L., Zhang, F., Stankiewicz, P., Halperin, J. J., Yang, C., Gehman, C., Guo, D., Irikat, R. K., Tom, W., Fantin, N. J., Muzny, D. M. & Gibbs, R. A. (2010) Wholegenome sequencing in a patient with Charcot-Marie-Tooth neuropathy. N Engl J Med 362, 1181–1191. Nakamura, K., Oshima, T., Morimoto, T., Ikeda, S., Yoshikawa, H., Shiwa, Y., Ishikawa, S., Linak, M. C., Hirai, A., Takahashi, H., Altaf-Ul-Amin, M., Ogasawara, N. & Kanaya, S. (2011) Sequence-specific error profile of Illumina sequencers. Nucleic Acids Res 39, e90. Nalls, M. A., Plagnol, V., Hernandez, D. G., Sharma, M., Sheerin, U. M., Saad, M., Simon-Sanchez, J., Schulte, C., Lesage, S., Sveinbjornsdottir, S., Stefansson, K., Martinez, M., Hardy, J., Heutink, P., Brice, A., Gasser, T., Singleton, A. B. & Wood, N. W. (2011) Imputation of sequence variants for identification of genetic risks for Parkinson’s disease: a meta-analysis of genome-wide association studies. Lancet 377, 641–649. Ng, S. B., Buckingham, K. J., Lee, C., Bigham, A. W., Tabor, H. K., Dent, K. M., Huff, C. D., Shannon, P. T., Jabs, E. W., Nickerson, D. A., Shendure, J. & Bamshad, M. J. (2010) Exome sequencing identifies the cause of a mendelian disorder. Nat Genet 42, 30–35. Nuytemans, K., Theuns, J., Cruts, M. & Van Broeckhoven, C. (2010) Genetic etiology of Parkinson disease associated with mutations in the SNCA, PARK2, PINK1, PARK7, and LRRK2 genes: a mutation update. Hum Mutat 31, 763–780. Pankratz, N., Beecham, G. W., Destefano, A. L., Dawson, T. M., Doheny, K. F., Factor, S. A., Hamza, T. H., Hung, A. Y., Hyman, B. T., Ivinson, A. J., Krainc, D., Latourelle, J. C., Clark, L. N., Marder, K., Martin, E. R., Mayeux, R., Ross, O. A., Scherzer, C. R., Simon, D. K., Tanner, C., Vance, J. M., Wszolek, Z. K., Zabetian, C. P., Myers, R. H., Payami, H., Scott, W. K. & Foroud, T. (2012) Meta-analysis of Parkinson’s disease: identification of a novel locus, RIT2. Ann Neurol 71, 370–384. Proukakis, C., Dudzik, C. G., Brier, T., Mackay, D. S., Cooper, J. M., Millhauser, G. L., Houlden, H. & Schapira, A. H. (2013) A novel alpha-synuclein missense mutation in Parkinson disease. Neurology 80, 1062–1064.

Annals of Human Genetics (2014) 78,243–252

251

L. Pihlstrøm et al.

Puschmann, A. (2013) Monogenic Parkinson’s disease and parkinsonism: clinical phenotypes and frequencies of known mutations. Parkinsonism Relat Disord 19, 407–415. Rivas, M. A., Beaudoin, M., Gardet, A., Stevens, C., Sharma, Y., Zhang, C. K., Boucher, G., Ripke, S., Ellinghaus, D., Burtt, N., Fennell, T., Kirby, A., Latiano, A., Goyette, P., Green, T., Halfvarson, J., Haritunians, T., Korn, J. M., Kuruvilla, F., Lagace, C., Neale, B., Lo, K. S., Schumm, P., Torkvist, L., Dubinsky, M. C., Brant, S. R., Silverberg, M. S., Duerr, R. H., Altshuler, D., Gabriel, S., Lettre, G., Franke, A., D’amato, M., Mcgovern, D. P., Cho, J. H., Rioux, J. D., Xavier, R. J. & Daly, M. J. (2011) Deep resequencing of GWAS loci identifies independent rare variants associated with inflammatory bowel disease. Nat Genet 43, 1066– 1073. Ross, O. A., Soto-Ortolaza, A. I., Heckman, M. G., Aasly, J. O., Abahuni, N., Annesi, G., Bacon, J. A., Bardien, S., Bozi, M., Brice, A., Brighina, L., Van Broeckhoven, C., Carr, J., ChartierHarlin, M. C., Dardiotis, E., Dickson, D. W., Diehl, N. N., Elbaz, A., Ferrarese, C., Ferraris, A., Fiske, B., Gibson, J. M., Gibson, R., Hadjigeorgiou, G. M., Hattori, N., Ioannidis, J. P., JasinskaMyga, B., Jeon, B. S., Kim, Y. J., Klein, C., Kruger, R., Kyratzi, E., Lesage, S., Lin, C. H., Lynch, T., Maraganore, D. M., Mellick, G. D., Mutez, E., Nilsson, C., Opala, G., Park, S. S., Puschmann, A., Quattrone, A., Sharma, M., Silburn, P. A., Sohn, Y. H., Stefanis, L., Tadic, V., Theuns, J., Tomiyama, H., Uitti, R. J., Valente, E. M., Van De Loo, S., Vassilatis, D. K., Vilarino-Guell, C., White, L. R., Wirdefeldt, K., Wszolek, Z. K., Wu, R. M. & Farrer, M. J. (2011) Association of LRRK2 exonic variants with susceptibility to Parkinson’s disease: a case-control study. Lancet Neurol 10, 898– 908. Ross, M. G., Russ, C., Costello, M., Hollinger, A., Lennon, N. J., Hegarty, R., Nusbaum, C. & Jaffe, D. B. (2013) Characterizing and measuring bias in sequence data. Genome Biol 14, R51. Ruark, E., Snape, K., Humburg, P., Loveday, C., Bajrami, I., Brough, R., Rodrigues, D. N., Renwick, A., Seal, S., Ramsay, E., Duarte Sdel, V., Rivas, M. A., Warren-Perry, M., Zachariou, A., Campion-Flora, A., Hanks, S., Murray, A., Ansari Pour, N., Douglas, J., Gregory, L., Rimmer, A., Walker, N. M., Yang, T. P., Adlard, J. W., Barwell, J., Berg, J., Brady, A. F., Brewer, C., Brice, G., Chapman, C., Cook, J., Davidson, R., Donaldson, A., Douglas, F., Eccles, D., Evans, D. G., Greenhalgh, L., Henderson, A., Izatt, L., Kumar, A., Lalloo, F., Miedzybrodzka, Z., Morrison, P. J., Paterson, J., Porteous, M., Rogers, M. T., Shanley, S., Walker, L., Gore, M., Houlston, R., Brown, M. A., Caufield, M. J., Deloukas, P., Mccarthy, M. I., Todd, J. A., Turnbull, C., Reis-Filho, J. S., Ashworth, A., Antoniou, A. C., Lord, C. J., Donnelly, P. & Rahman, N. (2013) Mosaic PPM1D mutations are associated with predisposition to breast and ovarian cancer. Nature 493, 406–410. Sauna, Z. E. & Kimchi-Sarfaty, C. (2011) Understanding the contribution of synonymous mutations to human disease. Nat Rev Genet 12, 683–691.

252

Annals of Human Genetics (2014) 78,243–252

Sikkema-Raddatz, B., Johansson, L. F., De Boer, E. N., Almomani, R., Boven, L. G., Van Den Berg, M. P., Van Spaendonck-Zwarts, K. Y., Van Tintelen, J. P., Sijmons, R. H., Jongbloed, J. D. & Sinke, R. J. (2013) Targeted next-generation sequencing can replace Sanger sequencing in clinical diagnostics. Hum Mutat 34, 1035– 1042. Vallania, F., Ramos, E., Cresci, S., Mitra, R. D. & Druley, T. E. (2012) Detection of rare genomic variants from pooled sequencing using SPLINTER. J Vis Exp 64, 3943. Vilarino-Guell, C., Wider, C., Ross, O. A., Dachsel, J. C., Kachergus, J. M., Lincoln, S. J., Soto-Ortolaza, A. I., Cobb, S. A., Wilhoite, G. J., Bacon, J. A., Behrouz, B., Melrose, H. L., Hentati, E., Puschmann, A., Evans, D. M., Conibear, E., Wasserman, W. W., Aasly, J. O., Burkhard, P.R., Djaldetti, R., Ghika, J., Hentati, F., Krygowska-Wajs, A., Lynch, T., Melamed, E., Rajput, A., Rajput, A. H., Solida, A., Wu, R. M., Uitti, R. J., Wszolek, Z. K., Vingerhoets, F. & Farrer, M. J. (2011) VPS35 mutations in Parkinson disease. Am J Hum Genet 89, 162– 167. Wang, K., Li, M. & Hakonarson, H. (2010) ANNOVAR: functional annotation of genetic variants from high-throughput sequencing data. Nucleic Acids Res 38, e164. Wang, S. R., Carmichael, H., Andrew, S. F., Miller, T. C., Moon, J. E., Derr, M. A., Hwa, V., Hirschhorn, J. N. & Dauber, A. (2013) Large-scale pooled next-generation sequencing of 1077 genes to identify genetic causes of short stature. J Clin Endocrinol Metab 98, E1428–E1437. Zimprich, A., Benet-Pages, A., Struhal, W., Graf, E., Eck, S. H., Offman, M. N., Haubenberger, D., Spielberger, S., Schulte, E. C., Lichtner, P., Rossle, S. C., Klopp, N., Wolf, E., Seppi, K., Pirker, W., Presslauer, S., Mollenhauer, B., Katzenschlager, R., Foki, T., Hotzy, C., Reinthaler, E., Harutyunyan, A., Kralovics, R., Peters, A., Zimprich, F., Brucke, T., Poewe, W., Auff, E., Trenkwalder, C., Rost, B., Ransmayr, G., Winkelmann, J., Meitinger, T. & Strom, T. M. (2011) A mutation in VPS35, encoding a subunit of the retromer complex, causes late-onset Parkinson disease. Am J Hum Genet 89, 168–175.

Supporting Information Additional Supporting Information may be found in the online version of this article at the publisher’s web site: Table S1 Bioinformatic pipeline. Table S2 Comparison between results from genotyping and pooled sequencing in 350 patients. Received: 13 December 2013 Accepted: 7 February 2014

 C

2014 John Wiley & Sons Ltd/University College London